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Preface 
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boundary  layer  flows.  The  study  was  accomplished  by  building 
upon  a computer  solution  method  developed  by  Dr.  Shang.  The 
study  allowed  me  to  apply  much  of  my  recent  academic  train- 
ing. 

I am  indebted  to  Dr.  Shang  and  Lt  John  Shea,  my  advisor, 
for  their  assistance  and  encouragement  in  all  phases  of  this 
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Abstract 

A numerical  solution  of  the  compressible  boundary 
layer  equations  was  developed  for  flows  over  either  two- 
dimensional  or  axisymmetric  surfaces.  The  solution  method 
is  an  extension  of  a computer  solution  developed  by  Dr.  J. 
S.  Shang  of  the  Air  Force  Flight  Dynamics  Laboratory.  The 
solution  method  is  capable  of  solving  for  boundary  layer 
parameters  in  either  laminar  or  turbulent  flows.  In  the 
case  of  turbulent  flow,  closure  is  achieved  by  use  of  a 
two-layered  eddy  viscosity  model.  The  boundary  layer  equa- 
tions are  solved  by  a numerical  marching  procedure.  A 
Mangler-Levy-Lees  transformation  of  independent  variables 
is  used  to  improve  the  efficiency  of  the  numerical  solu- 
tion. The  transformed  boundary  layer  equations  are  then 
linearized  by  a three  point  finite  difference  scheme.  The 
linearized  equations  are  solved  by  a matrix  solution  tech- 
nique. Comparisons  of  computed  boundary  layer  parameters 
with  experimentally  determined  parameters  were  made  for 
both  laminar  and  turbulent  flows  over  axisymmetric  bodies. 
The  comparisons  show  the  numerical  solution  to  be  very 
accurate. 
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NUMERICAL  SOLUTION  OF  THE  COMPRESSIBLE 


BOUNDARY  LAYER  EQUATIONS  OVER 


j 


AXISYMMEfRIC  SURFACES 


I . Introduction 


Purpose  of  the  Study 

The  purpose  of  this  study  was  to  develop  a numerical 
solution  for  compressible  boundary  layer  flow  over  axisym- 
metric  surfaces.  This  solution  was  to  be  applicable  to  both 
laminar  and  turbulent  flows.  This  study  was  accomplished  by 
expanding  an  existing  numerical  solution  procedure  developed 
by  Dr.  J.  S.  Shang  of  the  Air  Force  Flight  Dynamics  Laboratory. 

Existing  Numerical  Solution  Procedure.  The  existing 
numerical  solution  as  developed  by  Dr.  Shang  was  capable  of 
solving  for  boundary  layer  parameters  either  in  a two-dimen- 
sional case  or  the  limited  axisymmetric  case  of  conical  flow. 
This  numerical  solution  procedure  has  been  shown  to  yield 
accurate  results  for  compressible  boundary  layers  in  both  the 
laminar  and  turbulent  flow  cases  with  pressure  gradients  and 
heat  transfer  at  the  surface. 

The  original  algorithm  made  use  of  a transformation  of 
independent  variables  in  the  boundary  layer  equations  to  trans- 
form the  boundary  layer  equations  into  a coordinate  plane  in 
which  efficient  numerical  computation  was  possible.  This 
transformation  was  considered  only  in  two-dimensional  or 
conical  axisymmetric  flows  where  the  cone  radius  could  be  repre- 
sented  as  a linear  function  of  distance  along  the  surface. 

1 


The  transformation  also  employed  the  assumption  that  the 
boundary  layer  thickness  was  negligible  in  comparison  with 
the  radius  of  the  cone. 

The  transformed  boundary  layer  equations  were  linear- 
ized by  using  a finite  difference  method  of  approximating 
derivatives  within  the  equations.  These  equations  were  then 
solved  by  a matrix  method,  and  the  boundary  layer  parameters 
of  interest  were  obtained  by  an  inverse  transformation. 

Modified  Numerical  Solution  Procedure.  The  numerical 
solution  developed  in  this  report  follows  the  same  general 
method  of  solution  as  the  existing  solution  procedure.  The 
primary  difference  between  this  solution  and  the  existing 
solution  is  that  a more  general  transformation  of  indepen- 
dent variables  is  used.  By  generalizing  the  transformation, 
boundary  layer  flows  over  arbitrary  axisymmetric  bodies  can  be 
considered. 

The  transformation  of  independent  variables  used  in 
this  report  represents  the  axisymmetric  body  radius  as  a 
polynomial  function  of  distance  along  intervals  of  the  axis 
of  symmetry.  Using  this  model  of  the  radius,  the  transforma- 
tion of  independent  variables,  required  for  efficient  numer- 
ical computation,  can  be  applied  to  arbitrary  bodies  of 
revolution. 

The  new  solution  is  also  improved,  particularly  in 
the  case  of  slender  bodies  of  revolution,  by  accounting  for 
the  transverse  curvature  effects.  These  effects  were  ignored 
in  the  existing  solution  by  assuming  the  boundary  layer 
thickness  to  be  negligible  in  comparison  with  the  radius  of 


the  axisymmetric  body. 

The  boundary  layer  equations  were  then  transformed  by 
this  new  transformation,  and  the  resulting  equations  were 
linearized  by  using  a finite  difference  method  of  approxi- 
mating derivatives.  These  linearized  equations  were  then 
solved  by  the  same  matrix  method  of  solution  as  used  in  the 
existing  solution.  The  inverse  of  this  new  transformation 
produced  the  boundary  layer  parameters  of  interest. 

Outline  of  the  Report 

This  report  is  written  to  describe  the  modelling  and 
numerical  solution  of  the  boundary  layer  equations  over 
axisymmetric  or  two-dimensional  surfaces.  The  solution  pro- 
cedure described  in  this  report  was  implemented  in  a computer 
program  listed  in  Appendix  D. 

The  equations  necessary  to  model  the  boundary  layer 
flow  for  either  two-dimensional  or  axisymmetric  flows  are 
given  in  Chapter  II.  The  transformations  required  for  numer- 
ical computation  purposes  are  also  given  in  this  chapter. 

Chapter  III  provides  the  actual  mechanics  involved  in 
the  numerical  solution.  The  logic  of  this  numerical  solution 
as  implemented  in  the  computer  program  is  also  given. 

Chapter  IV  presents  the  results  of  two  test  cases  used 
to  validate  the  numerical  solution  procedure.  A comparison 
of  accuracy  with  two  other  numerical  solution  methods  is  also 
made. 

Chapter  V contains  conclusions  and  recommendations  re- 
sulting from  this  study. 
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II.  Mathematical  Formulation  of  the  Problem 


Governing  Equations 

The  boundary  layer  equations  describing  conservation 
of  mass,  momentum,  and  energy  as  developed  in  Appendix  A 
are 
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(2) 

(3) 


This  set  of  equations  can  be  used  to  describe  either  laminar 
or  turbulent  boundary  layers.  Tie  eddy  viscosity  terms  6 and 

A 

€ reduce  to  a value  of  unity  for  laminar  flow.  Also  the 

A*/  ■; 

normal  velocity  term  V,  which  is  the  sum  of  V and  , 

_ P 

reduces  to  v for  laminar  flow. 

The  superscript  j in  equations  (1),  (2),  and  (3) 
allows  the  equations  to  be  used  for  either  two-dimensional 
or  axisymmetric  flow.  For  the  two  dimensional  case,  o is 
set  equal  to  zero  to  remove  the  radius  r from  the  equa- 
tions. For  the  axisymmetric  case,  j is  given  the  value  of 
unity. 

_ A 

In  equations  (2)  and  (3),  € and  € are  considered  as 
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known  functions  of  P , u , and  y and  will  be  discussed 
later  in  this  chapter.  The  pressure  gradient  is  known  from 
the  inviscid  Euler  equation  (Ref  1^:21)  given  as 


ax  dx 


since  flow  properties  outside  the  boundary  layer  can  be 
determined  by  inviscid  flow  theory. 

Equations  (1),  (2),  and  (3)  then  contain  five  unknowns: 
q » P , if  • 4 t and  T . The  perfect  gas  law  and  Sutherland's 
viscosity  law  (Ref  17:19),  given  by  equations  (5)  and  (6) 
respectively,  are  used  to  express  //  and  P as  functions  of 
T and  thus  reduce  the  unknowns  to  and  T. 


P = P RT 


AL 


T«,W 
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The  pressure  Ms  a known  quantity  outside  the  boundary 
layer.  The  pressure  within  the  boundary  layer  can  be  consid- 
ered constant  in  a direction  normal  to  the  surface  bordering 
the  boundary  layer;  therefore,  equation  (5)  reduces  to  a 
relation  between  P and  T . 

The  use  of  equations  (4),  (5),  and  (6)  can,  therefore, 
be  seen  to  reduce  equations  (1),  (2),  and  (3)  to  a set  of 
three  equations  in  three  unknowns. 
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Transformation  of  Independent  Variables 


r 


1 


The  governing  equations  (1),  (2),  and  (3)  have  a singu- 
larity at  the  stagnation  point  ( r = O ) for  axisymmetric 
flow.  In  order  to  remove  the  singularity  and  also  to  remove 
most  of  the  variation  in  boundary  layer  thickness  along  the 
surface,  a Mangler-Levy-Lees  transformation  (Ref  2:261) 


/ \ 

?(X)  = 1%  U'AAt  ) 6X 

n(X,Y)=  (iiU  (r; 


[Ik  \ f lr)llUy 

(2 UJUj 


(7) 


(8) 


was  introduced.  The  removal  of  most  of  the  variation  in 
boundary  layer  thickness  improves  the  efficiency  of  numer- 
ical solution  of  the  boundary  layer  equations. 

The  quantities  r , rc , X , and  y in  equations  (7)  and 
(8)  are  those  shown  in  Fig.  1.  The  orthogonal  coordinate 
system  X,>  is  a surface  coordinate  system  in  which  the 


Fig.  1.  Surface  Coordinate  System 
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X axis  is  tangent  to  the  surface  and  points  in  the  stream- 
wise  direction  while  the  y axis  is  perpendicular  to  the 
surface.  The  term  r is  the  distance  from  the  centerline  to 
any  point  on  the  y axis  measured  perpendicular  to  the  axis 
of  symmetry.  The  term  lr  is  the  reference  length  used  in 
computing  the  free  stream  Reynolds  number.  The  ratio  £ in 
equation  (8)  is  known  as  the  transverse  curvature  parameter 
and  will  be  represented  as  t . 

The  partial  derivative  operators  associated  with 
equations  (7)  and  (8)  are 


Mt  Ue()[o_\  A.  \ -+-  ( i2ll  /jL  \ 
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The  governing  equations  (l),  (2),  and  (3)  are  now  transform- 
ed from  physical  surface  coordinates  X,y  to  the  transformed 
coordinates  '$,>1  . The  transformed  boundary  layer  equations 
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where  the  following  definitions  have  been  used. 
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When  equations  (5)  and  (6)  are  applied  to  equation  (16),  the 
expression  for  / becomes 
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Given  this  relation  for  / and  assuming  t » € . and  € can 
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be  expressed  as  functions  of  F , €> , and  \/ , equations  (11), 
(12),  and  (13)  become  a set  of  three  non-linear  differential 
equations  in  three  unknowns  F , © , and  V , The  expressions 

- A 

for  t , € , and  € are  developed  in  detail  later  irt  this 
chapter. 

Transverse  Curvature  Correction 

The  existing  numerical  solution  prior  to  modification 
assumed  that  the  transverse  curvature  parameter  t was  equal 
unity.  The  solution  developed  in  this  report  retains  the 
transverse  curvature  effect  in  the  solution  of  the  boundary 
layer  equations.  Since  the  parameter  t is  the  ratio  of  ^ , 
an  expression  for  t is  found  by  observing  from  Fig.  1 that 

r = -t  y cos  0 (21) 

Dividing  equation  (21)  by  £ and  using  the  definition  tr-  - 

'0 

gives 


t = • i + t 


(22) 


For  a given  value  of  X , the  derivative  of  t is 


dt  = 


Nr*) 


dX 


(23) 


Solving  equation  (23)  for  dX  and  substituting  the  value 
into  equation  (7)  yields 
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(24) 


When  equation  (24)  is  integrated  from  the  surface  outward 
in  the  normal  direction,  it  yields 


(25) 


Using  equations  (5)  and  (17),  the  expression  for  ~t  in  the 
transformed  coordinate  plane  is 


t 


2 U (Z ?)* ' 


CPS 
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(26) 


As  equation  (22)  indicates,  t = 1 is  a good  approxi- 
mation for  most  axisymmetric  boundary  layer  flows.  The 
numerical  solution  procedure  was,  therefore,  designed  such 
that  the  transverse  curvature  correction  could  be  retained 
or  dropped  from  the  solution  of  the  boundary  layer  equa- 
tions. 


Method  of  Modelling  Body  Radius 

The  geometry  of  an  axisymmetric  body  is  determined  in 

the  new  solution  procedure  by  fitting  a quadratic  curve  of 

« 

the  form 

f.(2)  = Ai  + B;Z  + 0,1*  (27) 


through  a series  of  surface  coordinates  as  shown  in  Fig.  2 . 
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Fig.  2.  Radius  Curve  Fitting 

Since  the  existing  computer  solution  procedure  required 
the  input  of  surface  pressure  at  a finite  number  of  surface 
locations,  the  coordinates  r0)Z  of  these  locations  are  input 
and  used  in  the  new  solution  procedure  to  determine  r;(2). 

The  coefficients  A;  , 8;,  and  C;  are  then  determined  by 
fitting  equation  (27)  through  three  consecutive  surface 
locations  starting  with  location  i . In  using  equation  (27) 
to  determine  f; , the  first  set  of  coefficients  A,.  , 6, , and 
C,  is  used  until  the  second  surface  location  is  reached  and 
then  the  second  set  of  coefficients  would  be  used  until  the 
third  surface  location  is  reached.  This  technique  is  contin- 
ued for  increasing  values  of  X until  the  last  surface  loca- 
tion ( i = A ) is  reached.  Since  the  changeover  in  coefficients 
is  made  at  points  for  which  r«  is  known  exactly,  errors  in 
determining  rs  do  not  propogate. 
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The  angle  0 between  the  centerline  axis  and  the 
tangent  to  the  surface  is  required  in  order  to  determine  the 
transverse  curvature  parameter  t in  equation  (26).  This 
angle  is  observed  from  Fig.  1 to  be 

0 - t^n"'  (jj-§)  (28) 

where  ^ determined  from  equation  (27)  is 


45  = B.*  + 

O t 


(29) 


Since  the  numerical  marching  technique  used  to  solve 
equations  (11),  (12),  and  (13)  proceeds  stepv/ise  along  the 
surface  of  the  body,  distance  along  the  surface  is  related 
to  centerline  distance  by  (See  Fig.  1) 


dZ  - Co  5 0 

d X 


Representation  of  Turbulence  Produced  Terms 

Closure  of  the  turbulence  correlations  in  the  boundary 
layer  equations  is  effected  by  the  eddy  viscosity  terms  € 
and  € (see  Appendix  A),  where 


e=l  + *- 

^ ‘ M 


e Pr 
+ «p~rt 


(32) 
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The  actual  eddy  viscosity  £ used  in  the  solution  procedure 
is  a two  layered  model  in  v/hich  two  separate  formulas  are 
used  for  the  inner  and  outer  region  of  the  boundary  layer. 

The  inner  region  model  is  based  on  Prandtl's  mixing 
length  theory  (Ref  13:  548). 


L\  = ?L1~  Jl 

* M 5 v 


This  expression  after  being  fitted  with  the  experimentally 
determined  mixing  length  4,  of  Van  Driest  (Ref  15:  1009)  is 

(£)']} |8|  (34] 

The  value  of  K,  used  in  this  report  is  0.4  as  given  by 
Van  Driest  and  validated  by  Cebeci  (Ref  4:  23).  The  wall 
shear  stress  4 is  given  as 


7w  ~ 


(gi 


The  eddy  viscosity  model  used  for  the  outer  region  of 
the  boundary  layer  (Ref  3 : 526)  is 

(S '!  = K*  •!  / K - w)  (36) 

o 

where  the  value  given  for  K2  is  0.0168. 

Changeover  from  the  inner  region  model  to  the  outer 


region  model  is  made  when  the  value  of  |-j.  is  equal  to 
Since  the  transformed  boundary  layer  equations  (11), 
(12),  and  (13)  are  solved  in  the  plane,  it  is  neces- 

sary to  transform  equations  (34)  and  (36)  into  this  coor- 
dinate plane.  These  equations  after  being  transformed  by 
use  of  the  transformation  equations  (9)  and  (10)  are 


P.  K, 
/ 


1 - exP 


-vnr 


24/ 


2 

an 

f (l  -f)  ©.  dh 
ea  Mi  l v t- 


(37) 

(38) 


where 

p, 

P 

The  turbulent  Prandtl  number  P^.t  appearing  in  equation 
(32)  is  determined  from  experimental  data,  and  the  value 
used  in  this  report  is  0.9  as  suggested  by  Cebeci  (Ref  2:260) 
for  air.  With  this  value  of  and  equations  (35)  and  (36), 

^ A 

£ and  € are  determined  since  the  static  Prandtl  number 
is  known  for  air. 

For  cases  of  transition  from  laminar  to  turbulent  flow, 
a transition  model  developed  by  Harris  (Ref  4:29)  is  used  in 
the  solution  procedure.  The  model  introduces  a transition 


r.t' 

a/'  e3 


= / 


W 


( 


an/v 


(39) 


(40) 
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factor  I which  is  multiplied  by  the  eddy  viscosity  € such 
that  equations  (3D  and  (32)  become 


(41) 

(42) 


The  transition  factor  I is  given  by 


l - exP 


-0.4-12  / y - y* 


(43) 


where  Xt  is  the  surface  location  at  which  the  transition 
from  laminar  to  turbulent  flow  begins.  With  this  model,  it 
can  be  seen  that  I would  have  a value  of  approximately  zero 
prior  to  the  transition  point  and  then  asymptotically 
approach  unity  beyond  the  transition  point. 


Boundary  Conditions 

The  transformed  boundary  layer  equations  (11),  (12), 
and  (13)  are  solved  by  a method  in  which  these  equations  are 
linearized  to  a first  order  set  of  equations  in  F , © , and 
v/  . This  set  of  equations  is  then  solved  subject  to  boundary 
conditions  on  F , © , and  \J  at  both  the  wall  and  the  outer 
edge  of  the  boundary  layer. 

The  boundary  conditions  at  the  wall  or  surface  of  the 
body  necessary  to  satisfy  the  conditions  of  no  slip  at  the 
wall,  no  mass  transfer,  and  a constant  body  temperature  are 

15 


- -- — ■ 


-a 


U ( Y,o)  = 0 

(44) 

v(x,o)  -0 

(45) 

T ( X,  o')  = 0 

(46) 

The  boundary  conditions  at  the 
layer,  which  are  determined  by 

outer  edge  of  the  boundary 
external  flow  properties,  are 

(4?) 

=Te(x) 

(43) 

In  terms  of  the  transformed  variables,  these  boundary 
conditions  are 

O 

ii 

o' 

U_ 

(49) 

v(/,o)  =o 

(50) 

e(?,o)  = t~/Tc 

(5D 

F (?,%}=  1 

(52) 

e(f.’u)  ='i 

(53) 

Although  there  is  no  physical  boundary  condition  for  V(X^), 
a mathematical  boundary  condition  on  V is  required 

in  the  numerical  solution.  Since  the  numerical  solution 


SS.TSPS 
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technique  uses  an  iterative  type  solution  at  each  step  along 
the  surface  to  obtain  a solution  to  the  transformed  boundary 
layer  equations,  the  value  used  for  a boundary  condition  on 
is  the  last  calculated  value  of  ^ (f , ) • 

- V(£,?u)i-.  <5*> 
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III.  Solution  Technique 

Grid  Generation  in  the  Transformed  Plane 

The  transformed  boundary  layer  equations  (11), (12),  and 
(13)  are  solved  by  a numerical  marching  procedure.  This 
procedure  obtains  a solution  at  each  step  along  the  surface 
moving  downstream  in  the  flow  field. 

In  order  to  implement  this  method  of  solution,  a grid 
is  generated  in  the  transformed  plane  in  such  a way  that  the 
ratio  of  any  two  successive  spacings  in  the  direction  is 
a constant.  Spacing  of  steps  in  the  ^ direction  is  variable 
in  the  sense  that  it  can  be  doubled  at  any  point  along  the 
surface.  The  grid  is  generated  by  the  following  formula. 

ATlr,.,  = (K)"'2  ATI,  n=2,N  (55) 


K is  the  ratio  of  any  two  successive  spacings  in  the 
direction  and  A^is  the  assigned  value  of  the  first  grid 
spacing  as  shown  in  Fig.  3*  This  grid  allows  for  short  steps 
near  the  wall  where  velocity  and  temperature  gradients  are 
normally  greatest  and  larger  steps  as  the  edge  of  the  bound- 
ary layer  is  approached.  The  values  of  K.  N.  and  A^are 
chosen  such  that  a minimum  number  of  points  have  to  be  used 
in  order  to  obtain  a solution  for  any  given  problem.  Typical 
values  of  K»  N,  and  Abused  in  this  study  were  1.1,  100, 
and  .0005  respectively.  A typical  starting  value  for  AX 
was  .001  which  was  then  doubled  at  several  points  along  the 
surface. 
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At  each  f station,  when  a solution  to  the  boundary 
layer  equations  is  obtained,  the  convergence  of  both  F and 
0 is  checked  at  the  (n  - 1 5)  grid  point  in  the  direction. 
This  is  accomplished  by  using  the  following  two  relations. 


F (M-IG)  - FCN-IS)  ^ O.OOOI  (56) 

©(M-16)  ~0(U-I5l  ~ 0.0001  (57) 

This  convergence  check  is  based  on  the  fact  that  f and  © 
approach  unity  as  the  boundary  layer  edge  is  approached.  If 
either  of  equations  (56)  or  (57)  is  true,  another  grid  point 
is  added  in  the  ^ direction.  This  solution  procedure  is 
then  capable  to  some  extent  to  correct  for  a value  of  N 
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A g \ F«hi  ,H‘i 
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9 

A 23  F***|  nti  0-vmt/ , K-l 

fl  ' K ' 

"F  022^  0 yn+i,  n 

023  0-KM+l.n+l  "'’0  21  + I K-l  C22  V-m+i.t! 

K ' *1  K 

02.3  Vo-w-n  K.+  I 

h 

' 

= D2 

•1 

(59) 

A 31  F on  + i K-l 

K 

A 32  F h 
n ' 

A33  Fw^(  rv+i  + 031  0»n+i.n-i 

r\  • n ' 

+ £>32  0*W+I,  * 

n 

+ 6>33n  0-yyx^t, , 

k + i + ©31  h-i  "^  032  V»v\+i  h 

n ' n 

033^  V»*t-n / r\ ti 

= 03 

h 

(60) 

The  subscripts  of  F , 6 , and  V indicate  the  grid  points  at 
which  the  variables  are  evaluated.  The  subscripted  coeffi- 
cients A, 6,C  and  D are  evaluated  from  known  values  of  the 
boundary  layer  properties  at  ^ stations  and  ^ as  shown 
in  Appendix  B. 


Matrix  Method  of  Solution 

The  linearized  boundary  layer  equations  given  in  equa- 
tions (58  )#(  59  )»  and  ( 60  ) were  developed  for  the  grid 
point  of  Fig. 3.  If  these  linearized  equations  are 

written  for  all  grid  points  of  £ station  1 with  M = 6, 

the  matrix  equation  (6l)  results.  This  equation  has  been 
written  for  an  unrealistic  case  of  six  grid  points  in  the 
direction  but  serves  to  illustrate  the  matrix  equation  that 
results  regardless  of  the  value  of  N . The  linearized  bound- 
ary layer  equations  have  fewer  terms  when  written  at  either 
grid  point  ,2  or  + n-1  . This  is  due  to  the  bound- 

ary conditions  which  are  used  to  specify  values  of  F , Q , 
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and  V at  grid  points  >n+i;l  and  + N . 

Matrix  equation  (61 ) is  then  a system  of  3x(N  -2) 
linear  equations  in  3x(N  -2)  unknowns.  Equation  (6l)  is 
solved  for  F , O , and  V by  use  of  a Gaussian  elimination 
method  (Ref  6:301). 

Once  the  solution  is  obtained,  it  is  checked  for  con- 
vergence by  computing 


F"  -m  -t 1 , 2 


*1  + 1,  I _ 


A*. 


= h 


COMV 


(62) 


each  time  equation  (61 ) is  solved.  The  present  value  and  the 
last  computed  value  of  FC0NV  are  retained  so  that  the  fol- 
lowing comparison  can  be  made. 


CO**V 


- F 

HCW 


CONV 


loco 


(63) 


The  value  of  ERROH  used  in  this  study  was  .0005.  If 
equation  (63 ) is  not  true  then  the  solution  of  equation  (61 ) 
is  accepted  as  the  solution  of  the  transformed  boundary 
layer  equations.  If  equation  (63)  is  true  then  the  sub- 
scripted coefficients  A , ft  * C , and  D are  updated  using 
the  values  of  F , © , and  V just  obtained  from  the  solution 
of  equation  (61  ) . With  these  new  values  of  A , 8 , C , and  D 
the  solution  of  equation  (61)  is  undertaken  again.  This 
procedure  of  updating  , B , C.  , and  D is  explained  in 
detail  in  Appendix  B,  This  iterative  type  procedure  is 
continued  until  the  convergence  criteria  of  equation  (63  ) 
is  satisfied  or  100  attempted  solutions  of  equation  (61  ) 
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Logic  of  the  Numerical  Solution  Procedure 

This  numerical  solution  procedure  as  implemented  into 
a computer  program  requires  certain  inputs  describing 
the  problem  in  order  to  obtain  a solution  to  the  boundary 
layer  equations.  The  free  stream  quantities  required  to  be 

input  to  the  computer  program  are  the  Mach  number  M , tern- 

00 

perature  , and  Reynolds  number  Re  . The  pressure  P at  a 
finite  number  of  surface  locations  and  the  surface  tempera- 
ture Tw .which  is  assumed  constant,  are  also  required  inputs 
to  the  computer  program.  It  is  also  necessary  to  know 
whether  the  flow  is  laminar  or  turbulent,  and  if  the  flow 
transitions  from  laminar  to  turbulent  flow,  the  transition 
point  X*  must  also  be  known.  These  quantities  necessary  to 
describe  the  physical  problem  are  depicted  in  Fig.  4. 


With  these  quantities  describing  the  problem,  the  curve 
fitting  procedure  to  determine  a continuous  function  for  re 
is  then  accomplished.  For  portions  of  the  surface  where  the 
radius  is  changing  rapidly,  the  number  of  surface  pressure 
entries  should  be  greater  than  on  portions  of  the  surface 
where  the  radius  is  slowly  varying.  This  improves  the  accu- 
racy of  the  radius  modelling  since  the  curve  fitting  proce- 
dure uses  the  coordinates  of  the  surface  pressure  entries  to 
determine  ro  as  a function  of  2 . 

In  order  to  make  the  transformation  of  independent 
variable  from  physical  coordinates  X,  y to  transformed 
coordinates  f it  is  also  necessary  to  know  the  local 
boundary  layer  edge  properties  of  and  Ue.  These  prop- 

erties are  determined  as  in  the  existing  solution  by  first 
determining  d?  at  each  of  the  surface  pressure  points.  The 
pressure  gradient  at  each  of  the  surface  pressure  points 
is  approximated  by  a central  finite  difference  method  as 
developed  in  Appendix  B.  This  approximation  for  ^ is 


dfM  = (X;  -X.-,)  P;»,  -X;-,  'X;*,')  P- 


~ (X;+,  -X;") 

' 7*J..  X;  -Xf-D. 


(64) 


where  the  subscripts  refer  to  the  surface  points  of 
Fig.  4.  The  pressure  gradient  was  then  made  a continuous 
function  of  distance  along  the  surface  by  a curve  fitting 
procedure  similar  to  that  used  for  the  radius  ^ .The  local 
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edge  pressure  can  then  be  determined  at  any  point  on  the 
surface  by  integrating  this  continuous  function  forward 
from  the  last  known  value  of  the  pressure  Pe.  . This  inte- 
gration is  accomplished  in  the  computer  program  by  use  of 
a trapezoidal  rule  for  integration.  The  boundary  layer  edge 
temperature  was  then  computed  by  assuming  that  temperature 
could  be  related  to  pressure  by  the  isentropic  relation 


(65) 


(Ref  9: 53)  where  and  Pe  are  stagnation  values  of  T and  P. 
With  these  quantities  known,  is  calculated  from  equation 
(6).  The  local  edge  velocity  is  determined  by  using  the 
energy  equation  for  an  inviscid  perfect  gas  (Ref  9:53)  which 
is  given  by 

CPT«  -Cp Te  (66) 

2 

Prom  the  perfect  gas  equation,  is  determined.  All  quan- 
tities required  for  the  transformation  equation  (7)  are  now 
known. 

The  first  step  AX  along  the  surface  is  now  taken  and 
transformed  into by  equation  (7).  The  grid  in  the  direc- 
tion is  also  generated  at  this  point. 

In  order  to  start  the  solution,  the  transformed  bound- 
ary layer  equations  are  solved  by  a similar  solution  tech- 
nique at  the  first  three  § stations. This  technique  assumes 
arbitrary,  but  identical,  profiles  for  F.0*  and  V at  the 
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first  two  5 stations.  Values  of  £ , € , and  t are  set  equal 
to  unity  for  this  similar  solution.  With  these  assumptions 
and  using  , all  the  coefficients  A , & , C , and  0 

of  the  linearized  boundary  layer  equations  ( 58  ) * ( 59),  and 
(60)  can  be  determined.  The  matrix  equation  (6l)  can  then 
be  solved  for  F , © , and  V/  at  £ station  3 of  Fig.  5*  When 
this  solution  is  obtained,  the  values  of  F , © , and  V at 
£ stations  1 and  2 are  set  equal  to  the  solved  values  at 
§ station  3.  The . coefficients  A , 0 , C , and  o can  now  be 
updated  and  a new  solution  for  F , 0,  and  \J  obtained  at 5* 
station  3.  This  procedure  is  repeated  until  the  solved 


values  of  F , © , and  V converge.  After  establishing  these 
initial  profiles  by  similar  solution,  the  solution  of  the 
boundary  layer  equations  can  proceed  as  described  in  the 


last  section.  Several  surface  steps  must  be  taken  before 
the  effect  of  starting  the  procedure  by  similar  solutions 
is  purged  from  the  solution  of  the  boundary  layer  equa- 
tions. From  comparison  with  experimental  data,  this  effect 
appears  to  have  vanished  within  the  first  twenty  steps 
along  the  surface. 

Once  a solution  to  the  transformed  boundary  layer 
equations  is  obtained  at  an  ^ station,  values  are  assigned 
to  F , 6 , and  V at  five  points  exterior  to  the  grid  used 
in  the  solution.  These  extra  values  of  F , 0 , and  V allow 
for  grid  expansion  in  the  'ty  direction,  if  necessary, 
later.  The  value  assigned  to  F and  6 at  these  extra 
points  is  unity.  The  value  assigned  to  V at  these  extra 
points  is  determined  from  an  expansion  of  V given  by 


V 


+ • j 


A 


N/-W  + I N -l 


(67) 


N,  N + 

where  can  be  determined  from  equation  (B-39),  and 

A^l^can  be  determined  from  equation  (55).  It  is  necessary 
to  calculate  values  of  F , 0,  and  V at  five  extra  points 
to  ensure  continuity  in  the  finite  difference  scheme  used 
to  approximate  derivatives  in  X • 

Having  solved  the  transformed  boundary  layer  equa- 
tions at  any  station,  the  boundary  layer  thickness  6 , 
displacement  thickness  6, , and  momentum  thickness  Sz  can 
be  determined.  The  boundary  layer  thickness  is  taken  as  the 
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value  of  ^ for  which  the  value  of  F is  approximately 
.995*  The  boundary  layer  thickness  in  physical  coordinates 
X , y is  determined  by  using  the  inverse  transformation  of 
equation  (8)  which  is 

y =XP-7 \ f -fv  (68) 

Z M'kt )J°  t 

The  integral  in  the  above  equation  is  evaluated  by  a trape- 
zoidal integration  from  the  surface  out  to  the  value  of  y\ 
for  which  F is  approximately  .995*  The  value  of  y deter- 
mined is  the  boundary  layer  thickness  & . The  displacement 
thickness  and  momentum  thickness  are  calculated  in  a simi- 
lar manner. 

_ A 

The  eddy  viscosity  terms  € and  £ are  now  calculated 
at  all  grid  points  of  the  ? station  for  which  the  solution 

of  the  transformed  boundary  layer  equations  was  determined. 

— A 

These  values  of  € and  € will  be  used  in  the  solution  of 
the  transformed  boundary  layer  equations  at  the  next  £ 
station. 

At  this  point,  f and  © are  checked  by  equations  (56) 
and  (57)  to  determine  if  they  are  approaching  the  outer 
edge  boundary  conditions.  If  necessary,  another  grid  point 
is  added  to  the  outer  edge  of  the  grid. 

Next  the  non-dimensional  heat  transfer  rate  and  skin 
friction  are  calculated  in  terms  of  a Stanton  number  and 
coefficient  of  friction.  The  Stanton  number  (Ref  16*526)  is 
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(69) 


St 


co 


He 

CP  Pm  <1* 


where  the  heat  transfer  coefficient  hc  is  given  by 


/A w Cr  /&T  A 

he  “ Pr(T.-UH^^ 


(70) 


The  coefficient  of  skin  friction  (Ref  13:128)  is 


Cf  = 

• "I o 

where  the  wall  shear  stress  7C*  is  given  by 


(71) 


(72) 


A local  Stanton  number  Ste  and  a local  coefficient  of 
skin  friction  Ct  are  also  calculated  using  equations  (69) 
and  (71)  with  P^  and  replaced  by  f*  and  respec- 
tively. The  wall  derivatives  j^^and  * appearing  in 

equations  (?o)  and  (72)  are  transformed  into  the  £ 
coordinate  plane  and  evaluated  by  a four  point  finite 
difference  scheme  developed  in  Appendix  C. 

At  this  point)  another  step  along  the  surface  could 
be  taken  and  the  procedure  for  solution  of  the  boundary 
layer  equations  repeated.  Several  other  parameters  are 
calculated  in  the  computer  program;  however,  the  logic 
presented  in  this  section  covers  the  main  points  of  the 
program.  A flow  diagram  is  given  in  Fig.  6 showing  the 
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logic  and  sequence  of  events  in  the  computer  program.  This 
flow  diagram  is  not  intended  to  be  a detailed  break-down  of 

i 

the  computer  program,  but  is  intended  to  show  the  sequence 
of  events  discussed  in  this  report. 


« 


Linearize  transformed  boundary 
layer  equations  with  finite 
difference  approximations 


Solve  linearized  equations  for 
• F , e , and  V at  all  grid  points 
of  the  present  f station 


Assign  values  to  F , e . and  V at 
five  points  exterior  to  grid 


Has  solution^---^^ 

of  the  linearizea'\^^^  No 

equations  converged  


Fig.  6.  Flow  Diagram  of  Computer  Program 


IV.  Validation  of  New  Solution 

Comparison  with  Experimental  Results 

The  two  cases  against  which  computed  solutions  v/ere 
compared  with  experimental  results  are  turbulent  flow  over 
a pointed,  waisted  body  of  revolution  and  laminar  flow  over 
a spherically  blunted  body  of  revolution.  These  two  cases 
test  the  new  solution  in  distinctly  different  types  of 
flows. 

Turbulent  Flow  Over  a Waisted  Body  of  Revolution.  The 
experimental  data,  upon  which  this  comparison  was  made,  was 
obtained  by  Winter  et  al.  (Ref  1 6 s 939 ) in  1965.  The  experi 
ment  provided  boundary  layer  data  for  turbulent  flow  over 
an  axisyrametric  body  in  the  presence  of  a streamwise  pres- 
sure gradient.  The  data  was  reduced  by  Winter  to  obtain 
velocity  profiles,  skin  friction  coefficient,  momentum 
thickness,  and  Mach  number  variation  at  the  boundary  layer 
edge.  The  accuracy  of  the  skin  friction  coefficients 
obtained  in  the  experiment  was  estimated  by  Winter  to  be 
about  t .02  per-cent.  This  was  the  only  mention  of  accuracy 
given  in  Ref.  16  for  the  experiment. 

The  geometry  of  the  test  body  of  revolution  is  shown 
in  Fig.  7.  Boundary  layer  data  was  experimentally  deter- 
mined for  several  different  free  stream  Mach  number  flows 
over  this  body.  A comparison  of  computed  boundary  layer 
parameters  with  experimentally  determined  determined  bound- 
ary layer  parameters  was  made  for  Mach  numbers  of  1 . 398  and 
1.70.  The  test  conditions  associated  with  these  two  Mach 
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numbers  were 


Case  1 Case  2 

M*  * 1.398  1.7 

R«  = 10.1  x 106  10.0  x 106 

40 

T«,  = 385. 388. 9°R 

~ = -976  .971 

* o 


Since  the  actual  pressure  data  obtained  by  Winter  et  al. 
was  not  available,  the  pressure  data  along  the  surface  of 
the  body  was  obtained  from  the  Mach  number  distribution 
determined  in  the  experiment.  This  Mach  number  distribution 
at  the  boundary  layer  edge  is  shown  in  Fig.  8.  The  pressure 
distribution  was  found  from  this  Mach  number  distribution 
by  using  the  following  perfect  gas  relation  (Ref  9 « 53 ) * 

Is  - ( I + *§■'  M*  j /V  ' (73) 

The  input  to  the  computer  program  required  the  pressure 
data  to  be  input  in  the  non-dimensional  form  ; there - 
fore  equation  (73)  was  multiplied  by  5*  to  obtain  the 
proper  pressure  input.  The  pressure  data  P»,  along  with 
corresponding  values  of  and  was  input  to  the 

computer  program  at  twenty-seven  surface  locations.  The 
reference  length  Lw  was  chosen  to  be  the  centerline 
length  of  the  axisymmetric  body. 

Comparisons  of  computed  and  experimentally  determined 
distributions  of  momentum  thickness  and  skin  .friction 
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mm. 


coefficient  are  given  in  figures  9 through  12.  Since  the 
axi symmetric  test  body  was  relatively  slender,  computed 
results  with  and  without  transverse  curvature  correction 
are  given  in  figures  9 through  12.  A comparison  of  computed 
and  experimentally  determined  velocity  profiles  is  given  in 
figures  13  through  15. 

These  comparisons  show  the  computed  boundary  layer 
parameters  to  agree  well  with  the  experimentally  determined 
parameters.  From  figures  9 through  12,  the  computed  solu- 
tion is  seen  to  be  improved  when  corrected  for  transverse 
curvature . 

Laminar  Flow  Over  a Spherically  Blunted  Cone.  The 
experimental  results,  upon  which  this  comparison  was  made, 
were  obtained  by  Bushnell  (Ref  1*18)  in  1968.  This  experi- 
ment provided  heat  transfer  data  for  laminarr  flow  over  a 
spherically  blunted  cone.  The  heat  transfer  data  obtained 
in  this  experiment  was  a ratio  of  heat  transfer  coeffi- 
cients where  h,  is  the  heat  transfer  coefficient  at 

ne  © 

the  stagnation  point.  The  heat  transfer  data  was  obtained 
by  thermocouples  and  an  automatic  data  reduction  system. 

The  estimated  accuracy  of  the  resultant  heat  transfer- 
coefficient  data  given  in  Ref  1 was  15  per-cent. 

The  blunted  cone  is  shown  in  Fig.  14.  The  test  condi- 
tions for  the  experiment  were 

= 7.95 

Re  = 1.65  x 106 

«0 

T*,  = 103.4  °R 
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The  pressure  data  obtained  in  the  experiment  was  given 

P* 

in  the  non-dimensional  form  . This  experimental  pressure 

» O 

data  was  not  given  for  the  spherically  blunted  portion  of 
the  cone;  therefore,  a Newtonian  pressure  distribution 
(Ref  10:366)  was  used  for  this  portion  of  the  test  body. 

This  combination  of  experimental  and  theoretical  pressure 
distribution  is  given  in  Fig.  17.  This  pressure  distribu- 
tion is  put  into  the  form  Rr  and  input  to  the  computer 

program  at  fourteen  surface  locations.  The  corresponding 
r 2. 

values  of  is  and  -r  are  also  input.  The  reference 
length  used  was  the  diameter  at  the  base  of  the  body. 

Since  the  heat  transfer  output  of  the  computer  program 
was  in  the  form  of  a Stanton  number,  it  could  be  compared 
with  the  experimental  data  by  use  of  the  following  equality: 


(74) 


The  stagnation  Stanton  number  was  in  error  due  to  the 
similar  solution  method  used  to  st$rt  the  numerical  solu- 
tion; therefore,  the  St  used  in  equation  (74)  was  found 

C 

by  a formula  given  by  Van  Driest  (Ref  10:366).  The  computed 

value  of  hs  determined  by  equation  (74)  is  compared  with 

hc-  u 

the  experimentally  determined  value  of  in  Fig.  18. 

The  computed  results  agree  very  well  with  the  experi- 
mental results  over  the  entire  surface  of  the  blunted  cone. 
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Comparison  of  Accuracy  l.Jith  Other  Computer  Solutions 


Since  reference  to  several  different  axisymmetric 
boundary  layer  programs  can  be  found  in  the  literature,  a 
comparison  was  made  with  two  of  these  computer  solutions 
to  determine  the  relative  accuracy  of  the  solution  develop- 
ed in  this  report.  The  two  computer  solutions  compared  were 
those  of  Cebeci  (Ref  2:370),  197^,  and  Harris  (Ref  k:  5^-55 )» 
1971.  Both  of  these  solution  methods  had  been  applied  to 
waist ed  body  data  of  Winter  et  al.  The  computed  momentum 
thickness  and  computed  skin  friction  coefficient  for  the 
^«r  1*398  case  were  used  as  the  basis  of  comparison. 

Cebeci' s method  of  solution  is  considerably  different 
from  the  solution  developed  in  this  report.  A different 
method  of  linearizing  the  boundary  layer  equations  is  used 
in  Cebeci' s method.  Also  in  this  method,  the  momentum  and 
energy  equations  are  solved  seperately,  and  an  iterative 
procedure  is  used  to  obtain  convergence  of  the  solutions 
of  these  two  equations. 

The  solution  method  of  Harris  differs  from  the  solu- 
tion developed  in  this  report  mainly  in  the  method  of  solv- 
ing the  boundary  layer  equations.  Harris's  method  solves 
the  momentum  and  energy  equations  simultaneously  and  uses 
these  results  to  integrate  the  continuity  equation. 

The  results  of  these  two  solution  methods  are  compared 
with  the  results  of  the  solution  developed  in  this  report 
in  Fig.  18  and  Fig.  19.  Since  Cebeci 's  solution  was  started 
using  an  experimentally  determined  velocity  profile  at 
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z/l.  = the  comparisons  with  both  solution  methods  was 
made  starting  at  this  point.  From  these  comparisons,  the 
solution  developed  in  this  report  appears  to  have  approxi- 
mately the  same  accuracy  as  Harris's  solution  and  slightly 
better  accuracy  than  Cebeci's  solution. 
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Momentum  Thickness  for  1.398 


O Experimental  Data 

— Numerical  Computation  with 

Transverse  Curvature  Correction 

Numerical  Computation  without 

Transverse  Curvature  Correction 


' .1  .2  .3  A .5  .6  .7  .8  .9  1.0 

Z/Lr 

Fig.  11.  Momentum  Thickness  for  M_=  1.7 

CD 

o Experimental  Data 

— Numerical  Computation  with 

Transverse  Curvature  Correction 

--  Numerical  Computation  without 
Transverse  Curvature  Correction 


Velocity 


Pressure  Distribution  on  Spherically  Blunted  Cone 


Dimensional  Heat  Transfer  Rate  on  Spherically  Blunted  Gone 


V.  Conclusions 


The  numerical  solution  of  boundary  layer  flows  as 
developed  by  Dr.  Shang  has  now  been  expanded  to  include 
flows  over  arbitrary  shaped  axisymmetric  bodies.  The  new 
solution  was  validated  by  testing  it  against  two  different 
sets  of  experimental  data.  The  experimental  data  of  Winter 
et  al.  (Ref  16:939)  tested  the  computed  values  of  momentum 
thickness,  skin  friction  coefficient,  and  velocity  profiles 
for  turbulent  flow  over  a slender  pointed  body  of  revolu- 
tion. The  data  of  Bushnell  (Ref  1:18)  was  used  to  test  the 
computed  heat  transfer  rate  for  laminar  flow  over  a spheri- 
cally blunted  cone.  The  comparison  of  these  computed  bound- 
ary layer  parameters  with  the  experimentally  determined 
parameters  show  the  new  solution  to  be  very  accurate. 

The  accuracy  of  the  solution  procedure  developed  in 
this  report  was  compared  with  two  other  recent  computer 
solutions  using  the  data  of  Winter  et  al.  The  solution 
developed  in  this  report  was  determined  to  have  approxi- 
mately the  same  accuracy  as  Harris's  method  of  solution 
(Ref  4:5^)  and  slightly  better  accuracy  than  Cebeci's 
method  (Ref  2:370). 

The  solution  procedure  developed  in  this  report  can, 
therefore,  be  used  to  accurately  determine  boundary  layer 
parameters  for  either  laminar  or  turbulent  flow  over  either 
a two-dimensional  or  axisymmetric  surface. 
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Appendix  A 


Development  of  Two-Dimensional  and  Axisvmmetric  Boundary 

Layer  Equations 

The  boundary  layer  equations  are  obtained  from  equa- 
tions describing  the  conservation  of  mass,  momentum,  and 
energy  within  a moving  fluid.  These  conservation  equations 
can  be  obtained  by  examining  the  rate  of  change  of  mass, 
momentum,  and  energy  within  a finite  region  R of  space 
enclosed  by  a surface  S fixed  in  the  flowing  fluid.  Using 

A 

Fig,  21  in  which  n is  a unit  vector  normal  to  the  surface 
5 . the  integral  equations  for  time  rate  of  change  of 
mass,  momentum,  and  energy  (Ref  11:26-27)  are  respectively 


Fig,  21.  Finite  Region  of  Flow  Field 
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TtJJf  Pit  = -ffpv  *n  d5 

sjT^Vdr  = -ffyipy  *n)dS 

+^g-h  ds 

Jtj[  P(et^)dr»'/(°(e+^)v-ndS 

R 

4-JT  (g-  -vVhdS  ~<f  A VT*n  d5 


(A-l) 


(A-2) 


(A-3) 


Equations  (A-l),  (A-2),  and(A-3)  are  valid  for  compress- 
ible, viscous,  heat  conducting  fluids.  The  vector  is 

the  heat  flux  vector  and  g*  is  the  stress  tensor  (Ref  11:35) 
given  by 


O'  = -PI,  -V  (XV*  v)Td  +^(VV+VV*)  (A-4) 

The  value  of  A in  equation  (A-4)  as  determined  by  Stokes 
(Ref  13:57)  is  -f*. 

By  use  of  the  divergence  theorem  (Ref  5:292)  as  given 
in  equation  (A-5),  the  surface  integrals  in  equations  (A-l), 
(A-2),  and  (A— 3 ) can  be  changed  to  volume  integrals. 


ff  A • fi  dS  =fff  V • AdT 


(A-5) 


Since  R.  is  a fixed  volume,  the  partial  derivative  can 
be  taken  inside  the  integral  sign  of  equations  (A-l),  (A-2), 
and  (A-3).  By  combining  terms,  each  of  these  equations  could 
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be  expressed  as  one  volume  integral.  By  taking  the  limit 
as  R tends  to  zero,  the  differential  form  of  these  conser- 
vation equations  is 

+ KJ-PV  - 0 (A-6) 

dt  * — 

) + V ■ (f  V V -£  ) = ° (A-7) 

+ V-AT  + -V)  (A-8) 

In  order  to  further  simplify  equations  (A-6),  (A-7),  and 
(A-8),  the  assumption  of  steady  flow  is  made.  Also  the 
equation  for  mechanical  energy,  obtained  by  forming  the 
scalar  product  of  equation  (A-7)  with  the  velocity  vector 
V , is  subtracted  from  equation  (A-8).  With  these  simpli- 
fications, equations  (A-6),  (A-7),  and  (A-8)  become 


V • PV  =0 

(A-9) 

v-(  i°yv  -q)  = o 

( A-10) 

V • fhV  ^V-VP  tV'^VT  +•  f -vy. 

(A-ll) 

The  definitions  of  enthalpy  ( h = e t ^ ) and  rate  of  strain 
tensor  ( jf  = -PI^  ) have  been  used  in  equation  (A-ll). 

In  order  to  obtain  boundary  layer  equations  from  equa- 
tions (A-9),  (A-10),  and  (A-ll),  it  is  necessary  to  express 
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these  equations  in  a surface  oriented  coordinate  system. 

The  coordinate  system  to  be  used  for  an  axisymmetric  sur- 
face is  shown  in  Fig.  22  .This  orthogonal  curvilinear  coor- 
dinate system  is  oriented  such  that,  for  any  point  on  the 
surface,  the  X axis  is  parallel  to  the  surface  and  the  y 
axis  is  perpendicular  to  the  surface.  The  third  orthogonal 
coordinate  uj  is  taken  as  the  azimuth  angle  within  the 
meridian  plane  of  the  axisymmetric  body.  The  quantity  Kc 
is  the  curvature  in  the  £,?  plane  ?and  is  the  radius 
of  curvature  within  this  plane.  The  curvature  Kc  is  defined 
to  be  ^ . The  scale  factors  (Ref  12 <220)  necessary  to 


Fig.  22.  Orthogonal  Curvilinear  Coordinates  for  an 
Axisymmetric  Body 
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define  differential  arc  lengths  along  any  one  of  the  X,  Y, 
or  w surfaces  while  holding  the  other  two  coordinates  con- 
stant are  given  by  h„,  >V  and  hw  of  the  following  three 
relations t 


Mx  = ( 1 + Kcy ) jx 

(A-12) 

Mr  My 

(A-13) 

Kudto  = rjw 

(A-14) 

With  these  values  of  hx=(+Kcy,  hY=i  , and  V , the 
vector  operators  in  the  surface  oriented  coordinate  system 
can  be  obtained  by  using  the  general  form  of  these  opera- 
tors (Ref  5*302-303)  which  are 


V*A  “J. 


h.h,  h3 


i.  (h.hjA,) (h,hjAa)+  (h,h2A  jl  (A-15) 
3X(  aXz  8X3  J 


r-7  U/  _ I djf  ^ 4.  I ^ ' \ l?  2 

~ h,  dXi  hi  h%  <)Xj  * 


(A-16) 


V X A = I 


h,  he  h3 


ht  ? 

iy, 


A 

j 

* 

TX2 


K.Ai  htAj 


A 

H,  h 


\x3 


h$A3 


(A-17) 


In  the  above  three  equations,  X, , X2,  and  X3  are  the  orthog- 
onal curvilinear  coordinates  along  which  the  unit  vectors 

A A A 

i » 0 » and  A are  respectively  aligned. 


Before  using  equations  (A-15)*  (A-l6),  and  (A-17)  to 
form  vector  operator  for  the  orthogonal  curvilinear  coor- 
dinates of  Fig.  22,  the  following  approximation  is  made. 

| + Kc  Y ^ 1 (A-18) 

This  is  a valid  approximation  since  the  maximum  value  of  y 
will  be  the  boundary  layer  thickness  and  the  magnitude  of 
Kt  is  small  everywhere -on  an  &xisymmetric  body  except  at  a 
sharply  pointed  tip.  However,  the  value  of  y is  approxi- 
mately zero  at  the  tip  since  the  boundary  layer  thickness 
is  very  small  in  this  region.  With  this  assumption,  the 
vector  operators  for  the  orthogonal  curvilinear  coordinates 
of  Fig.  22  are 

V-A  = J_fi  (rAxUi  (rAv)  +i_  (Ajl  (a-19) 
r L d*  »y  jSi  J 

Vf  = 13?  * * j + (a-2o) 

ax  ay  1 d 

y x a = j_ 

r 

These  vector  operators  are  used  in  equations  (A-9)*  (A-10), 
and  (A-ll)  to  express  the  conservation  equations  in  a sur- 
face oriented  coordinate  system. 


i 

Jl 

ax 


j. 

av 

Av 


rk 

a 

Tz 


(A-21 ) 
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Boundary  layer  equations  applicable  to  turbulent  flows 
over  axisymmetric  surfaces  can  be  obtained  by  the  following 
steps i 

1.  Apply  the  vector  operators  of  equations  (A-19)> 
(A-20),  and  (A-21)  to  the  conservation  equations 
(A-9) i ( A-10) , and  (A-ll) 

2.  Assume  no  swirl  effect  ( -§w  = o ) 

3.  Represent  the  fluid  properties  as  the  sum  of  a 
mean  flow  property  and  a fluctuating  property 
(Ref  13*526) 

4.  Time  average  the  resulting  equations  and  retain 
only  the  significant  correlations  (Ref  4:13) 
which  are  P'xr',  aV|and  'U'T' 

5.  Perform  an  order  of  magnitude  analysis  using  the 
thin  boundary  layer  assumption  (Ref  13*118). 

The  boundary  layer  equations  for  conservation  of  mass, 
momentum,  and  energy  which  result  from  this  procedure  are 
respectively 

Js,[r?u]  + l»  \_f  P (v  + e-f  )]  r o (A-22) 


(“«  lx  p ) Sy  - ‘ Sx 


+ r iy[r(# 


(A-23) 


I 


t^')£  = « ft 

-fW'  ft  + Uy(^  ft)  ' ft 

~f  i (^  ft)  ' r iv{rcF?  iFF) 


(A-24) 


The  radius  entered  into  equations  (A-22),  (A-23)» 
and  (A -24)  by  the  scale  factor  r . For  a two-dimension- 
al surface,  the  values  of  W*  and  hy  would  be  unity,  and  the 
value  of  hw  would  be  zero.  The  result  of  deriving  boundary 
layer  equations  for  the  two-dimensional  case  would  be  equa- 
tions (A-22),  (A-23).  and  (A-24)  with  T replaced  by  unity. 
In  order  to  use  the  same  set  of  equations  for  either  two- 
dimensional  or  axisymmetric  surfaces,  P is  replaced  by  TJ 
where  j =0  for  two-dimensional  surfaces  and  i = I for 
axisymmetric  surfaces. 

In  order  to  have  a solvable  set  of  boundary  layer 
equations,  ^ is  grouped  together  with  V,  and  the  turbu- 
lent transport  quantities  P U‘V and  V"['  are  represent- 

ed by  an  eddy  viscosity  and  eddy  conductivity  respectively. 
This  method  of  handling  the  turbulence  produced  correla- 
tions is  introduced  into  equations  (A-22),  (A-23),  and 
(A-24)  by  use  of  the  following  definitions: 


— P’V 

v = v + —r 


(Normal  velocity) 


(A-25) 


(Eddy  viscosity) 


(A-26) 
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Kt  = 

"Co  ? WV 

35 

Ty 

(Eddy  conductivity) 

(A-27) 

Prt  = 

Cp6 

K« 

(Turbulent  Prandtl 
number) 

(A-28) 

€ = 

1 + i ) 

(A-29) 

A 

€ = 

1 

1 + % 

t 

(Associated  eddy 
’ viscosity  terms) 

(A-30) 

After  replacing  r by  rJ  and  using  the  above  defin- 
tions,  equations  (A-22),  (A-23),  and  (A -24)  are 


fx  (r^u)  + ly  (rV#)  = 0 u-31) 

+pv  57=-|-y  + rJ *Y {^as  ly)  (A-32> 

PU7*(c^)*isV&fcfi)=“  $ + *?($)* 

+ Mr(^*v(c,T))  u.,„ 

This  set  of  boundary  layer  equations  (A-31),  (A-32),  and 
( A— 33 ) is  valid  for  compressible  heat  conducting  boundary 
layers  over  either  two-dimensional  or  axisymmetric  surfaces 
depending  upon  the  choice  of  j . These  equations  are  also 
valid  for  either  laminar  or  turbulent  flows  depending  upon 

_ A 

the  value  given  € and  € . 


Appendix  B 


Three  Point  Finite  Difference  Approximations  and  Li ear- 
ization  of  the  Transformed  Boundary  Layer  liquations 


The  foundation  of  finite  difference  methods  is  the 
Taylor  series.  Both  the  forward  difference  method  used  to 
approximate  derivatives  in  § and  the  central  difference 
method  used  to  approximate  derivatives  in  are  obtained 
by  manipulation  of  a Taylor  series.  These  approximations 
for  derivatives  are  used  to  linearize  the  transformed 
boundary  layer  equations.  The  finite  difference  approxi- 
mations are  to  be  determined  for  a grid  point  'yn*\iv\  of 
Fig.  3. 


The  forward  difference  method  is  obtained  from  the 
following  two  truncated  Taylor  series  expansions  of  a 
variable  H about  the  grid  point  »+i,n  . 


- H 


n 


-AS. 


id 

as 


2 W 


(B-l ) 


2 as* 


( B-2) 


*»*♦»,  h 


Solving  equation  (B-l)  for  and  substituting  the  value 

into  equation  (B-2)  yields 
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'a?,  + ttu  ' 


H 


| Ag.  HCi 

l Af.At, 


+■  f"DCa  ri*-i(  h 


(B-3) 


The  central  difference  method,  used  to  approximate 
derivatives  in  V(i  is  obtained  from  the  following  two  trun- 
cated Taylor  series  expansions  about  a grid  point  >»♦' , 


H *v»i,  <*■‘1 


H V|4l(n-I 


Hx+i,n 


H -Vl  + >,  M 


-4-  m 

A** 

4-  l 

(B-4) 

an 

z 

a^| 

- c)H 

1 i.., 

+ i. 

£Ua 

an! 

(B-5) 

4h4>lfM 

z 

i*U 

Solving  equation  (B-4)  for  and  substituting  the  value 

into  equation  (B-3)  gives 


an 

U>ln-. 

1 

H 

'A>Ih-.  -/Pi* 

a»v 

j ^-Hl+1,  H 

- 1 

[A?U 

t — \ 

1 H Kl't 

( 

[A»V,(A?U-.+A*Oj 


Since  the  transformed  boundary  layer  equations  contain 

d*H 

second  derivatives  in  H , the  approximation  for  is 

found  by  eliminating  ib  between  equations  (B-4)  and 

dh. 

(B-5). 
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1m  -\p  In 

<>  n*  |A>ln(A**  + A>U-.)J  ATI 


+ 1 


(B-7) 


In  order  to  linearize  products  of  variables 
it  is  first  necessary  to  find  expressions  for  both  H and 
at  G-  station  -m  + j as  functions  of  the  known  values  at  5 
stations  and  Vi  . The  desired  expression  for  H,*,*,,*  is 

found  by  eliminating  ^ between  equations  (B-l)  and  (B-2). 
After  truncating  second  order  terms,  the  expression  for 
is 


H^,,n  = *U,„  - [Mil  H 

L j lA*,J 


(B-8) 


Replacing  H by  G-  in  equation  (B-8)  gives  the  desired 
expression  for  G . When  equation  (B-8)  is  used  to  eval- 
uate and  , the  values  will  be  given  by  mm  and 

Gmi  respectively.  The  following  equality  is  now  used  in  the 
linearization  process. 


H Gt-hi+  (jn  — + H+w-tijti  GVwfi,  h “ i (ntrrtij  « (g— 9) 

By  replacing  terms  in  equation  (B-9)  with  values  determined 
from  equation  (B-9),  the  following  expression  linear  in 
and  G-  is  obtained. 

If  n 
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n O’ - 


(B-10) 


HMI  (?w+((h  - HMI 


1 (r  M I 


Products  of  variables  and  derivatives  (G-tJL,  or  G-  f£|  ) 

can  now  be  put  into  linear  form.  Using  equation  (B-3),  the 
product  G is 


= ^ h C 

H h 

A?, 

[as. 

Of,  6?, 

H, 


+ 


4£z 


(B-ll ) 


Equation  (B-10)  is  used  to  replace  the  product  h „ 
in  equation  (B-ll),  which  makes  the  equation  linear  in 
and  . if  the  evaluation  of  ^ | ^ is  given  by  HV, 

the  linear  approximation  for  G is 


dn 


= G 


'HV  + 


G 


- GMi  • HY 


(B-12) 


A product  of  derivatives  in  n is  linearized 


by 


IW) 


*»l*l,* 


r 


hV-(^L,„  -SY-Hy  (B-13) 


where  &y  is  the  evaluation  of  ^ 

_ . ■»»,  ti 

Before  using  these  approximations  for  derivatives  and 
product  terms  to  linearize  the  transformed  boundary  layer 
equations,  the  following  definitions  are  introduced  to 
simplify  the  finite  difference  expressions. 
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11 


v,  = a Ur,  4-  a 
Ar,  + btz 

Yz  - 2 (flg,  +afj 

a*. 

Xjr  a ar/ 


a?,  u?,  + a?£) 

Y* = A ? . +•  a f? 

Ar, 

Y«--  AJ, 

Af, 

Y7  - A^h 
Y«=  2^! 


a>iM-.  (^v^CT) 


Y*  = 2 A>u„ 


+ an, 


Y,o  “ A^i -ahm 
AAK_, 


(B-14) 

(B-15) 

(B-16) 

( B— 1 7 ) 
(B-18) 

(B-19) 

(B-20) 

(B-21) 

(B-22) 


Ml 


tv  - ^ 

• v - a* 


(B-29) 


VMl  - YfV*,,*  - Ys-  V*,-.,  M 


(B-30) 


VM2  ~ /sV-**,*  *”  V3  V****- i/M 


(B-31) 


(B-32) 


The  transformed  boundary  layer  equations  (11),  (12),  and  (13) 


with  derivatives  expanded  are 


F +-  ■?  r ^ ^ = 0 

r f 25  at  In  u 


z_  .*i  A 3?  a( t2J)  ii 

/ £ t ^ + / & ax  an*  + ft  an  an 


(B-33) 


+ et*‘j£&  & -vft  ^^e-(3FJ  *2JF$  (B-3M 

/ I t.  £ bCt  ) d&  + JL  &£.  t &© 

pP  an*  Pr  an  an  pr  >n  an 

+/  4/  fd©f£t2J  - Vi®  + /€  t*J«*  (f£) 

p„  a©  vanj  an 


= 2?F^ 


(B-35) 


r 


Using  the  development  for  linearizing  terms  in  the  variables 
Cj-  and  H , the  transformed  boundary  layer  equations  can  now 
be  linearized  in  F ,e,  and  V . At  the  grid  point  + , 

the  linearization  is  accomplished  by  substituting  the  follow- 
ing expressions  into  equations  (B-33),  (B-34),  and  (B-35). 


af 

r 

Y|  F n 

Y?  F->vi,  n + Y5  F 

(B-36) 

a? 

2 

6f 

“Vi  ("wH,  n-l 

" Y10  + + Y9  Emti,  nti 

(B-37) 

46 

_ 

~Yg 

r 

Y*©  Y9  nf  1 

(B-38) 

a* 

'"1+1,  * 

2 

ATi* 

&L 

Yg  \f/ry\+i.  <*-l 

Yio  Y*v,+t  ^ + Y?  V/M-t+r,  »>•(•) 

(B-39) 

w 

/>h+i,k 

a*F 

av 

— 

Yg 

1 ^ ^7  F*n-*->,K  "+  Yt  F’/mH-l/K*l 

(B-40) 

&K 

a** 

_ 

Yfc  H- 

1 ^ Y7  Q'tyt-ri.  +”  Y$  Q-m+l  *i+< 

(B-41) 

a*; 

fifVae'il 

- FY/dS^  + TY  /4J\  - FY-GY 

,k  \^J*n-n,h  [d?! 

(B-42) 

89 

2 

1*1 

= 2 FY(^ 

- ^y2 

fc, 

(B-43) 

m 

f 

_ a ty  ^ 

A - TY2 

(B-44) 

\a>\, 

F«|-  » 

FM2)F„,,„  -Fnia-y, 

(B-45) 

* 


I -m-f  i,r» 


a or. 
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'! 

, 


Fae 

— 

(y,  * TM!  -TM2)  F^+,,h  + Y,  *©^4)jK 

h 

2Af2 

-Y,  • TMI  - FM\ 
ZiSK-t 

(B-46) 

- VMI-FY 

V 9'?1//m+  iy  h 

(B-47) 

V 

&y\ 

"W+  l,h 

14*. +.,*  'TV  + -VHI'TV 

N /*M  + i,H 

(B-48) 

F2 1 

- 

2 * FMi-  F^.k  - F™l* 

(B-^9) 

w+ 1,*, 


After  collecting  terms,  the  result  of  substituting  equations 
(B-3 6)  through  (B-^9)  into  the  transformed  boundary  layer 
equations  (B-33).  (B-3^),  and  (B-35)  is  the  following  three 
equations. 

F^-ri,*  +•  \4n  + M-,  £ - Ys  ' XL.J 

+ -Xl]  = £ -FM2(B-50) 

i 


'/0 


*m+\ , *•' 


|Y>  XL  [2  /€t~  - /£  it 

+ vm\  - € t20  il  Tyj 

de  'J 

F**.,*  ["  ^VXL/i  t2J  -2  Y.c-XL  (/€ 

+ /t2J£f  -vmi  + e t2J-M  -TY  " 2^« '^' 
' Zb%z  - (3  f HI  (.2  Y,  • FHI-FM2') 


l'»n-H<  1 


4 


T .20 


2Yi0  ' XL • £ t 


*/ 

de 


FY 


4-01H+1,  i-i 


+ 0-»«+i<i+« 


■f  V/m+i,  i 


[-Y.-XL-if^  FV 

[ y?  'XL  * £ t2Jai  fx 

L a©  J 

-FY' Aft]  = [fy-AJ2(6  t2J 


^ TY 
de 


-FW\2(A  + vO 


vwi) 

(B-5 


,„.k-,[-2  VFY-XL-/e  tai' 


v,0  FY  • XL'  / € t!J“  -fdV.'TM|-TM2l 


+ F^,,*„  I 2 Y,  -FV  XL'/e  t"* 

+ 0->h+ij^-i  Vb  XL  (2  / fc  t -/€  - 

. pr  \ 


. 2 J 

d€  t 


tpr  • VMl  -2Ty  U i t2J 

/ 


y/b\+\,h 


-Y  Y7  “XL  /6t2J  - £ Yco  ( xl') 

pr 


f/fi 
\ 


+ /-il  -t23  -Pr*VMl  + 2TY  i/  e tai| 


, xl  Cav4  /e  * 

' L Pr  t 


+ v,  f/€  + / ai  t2J 

N ^1  a?i 


-fV  * VMl  + 2 TY  4 t2J 

a©  ' 


+ v~,+1,  k ty]  ■=  fl?g-ry  lue  a t”  ry 

- Pp  \ 

-Pr  WW&?,  /?  t*5  •‘•FYa-£'Y,'TMl-FI*U  ( B- 52 ) 

The  derivatives  in  7l  , contained  in  the  brackets  []  ] of 

equations  (B-50),  (B-51),  and  (B-52),  are  evaluated  using 
knovm  values  at  % station  7* . The  derivative  is  found 
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by  explicitly  taking  the  derivative  of  equation  (20)  and 
evaluating  the  result  at  f station  >*. 

ti  ^ a 

The  values  of  t , f , e . and  e within  the  brackets  are 
determined  from  known  values  at  f stations  •»»-•  and  '*<  by  use 
of  the  relation  given  in  equation  (B-8).  All  other  quantities 
in  the  brackets  can  be  explicitly  determined  at  the  grid 
point  n . Equations  (B-50),  (B-51).  and  (B-52)  are 

then,  at  a given  f station,  a constant  coefficient  set  of 
linear  algebraic  equations. 

The  solution  method  for  these  equations  uses  an  itera- 
tive method,  if  required,  to  achieve  the  convergence  criteria 
of  equation  (63).  If  iterations  are  required,  the  values  of 
F Ml  » TM l » VK|  , FV  . TV  , and  vy  are  updated  after  the 
initial  solution  of  equations  (B-50),  (B-51).  and  (B-52)  at  a 
given  % station.  The  values  of  FMi  , TMl  , and  VM|  are 
replaced  by  the  values  of  F , 6 , and  V determined  in  the 
previous  iteration.  With  these  values  of  F , 6>  , and  V,  the 
derivatives  ^ and  can  be  determined.  FV,TV» 

and  vy  are  then  replaced  by  If  , , and  . 

Equations  (B-50),  (B-51)*  and  (B-52)  are  now  written 
in  a simplified  form  for  reference  purposes.  Equations  (B-50), 
(B-51).  and  (B-52)  are  given  by  equations  (B-53).  (B-54),  and 
(B-55)  respectively. 
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V Fw  + i,*-!  F/»»ft  i,  *1  + An  F*m4  J,  h + l t"  ft  (,  0 -*»»•»•  i 7 ki-/ 


+ 6l2  0**+,,*  +613  « t-l  C||  ft-l  *t  C'lZ  1 

n n ' n In  ' 


Cl3  V-mtl  n+i  D 1 

* » n 


(B-53) 


A 2 1 F»,+i,  n-i  t-  A 22  _ F 1,  h + A23  F**i-t  I ItH  Q'y*  +1  n - ) 

n » n 1 ' n 7 


t"  &22h  0-**»+i , 11  6 33„  ©'**»■*■(, h + i "t  ^~2l M V>,+|  _ n-i  ^2?  \/»*\+  1, 


V-M+ijK+i  ~ 0; 


K 


(b-54) 


A31  F*,+  1;  n_,  + A 3 ? F-**-t  i,  n Ajj  F^,+l>  ■*■  C*3iM  ©'m-»i/n-i 


+ &J2  O'*!*!,  n + ®33  + ^3,h  Vtrn-,'H-i  ''"  ^-32k  V^t  + i, 


^"33^  V^,4i,  uti  _ 


(B-55) 


The  subscripted  coefficients  A . 8 » C , and  0 are  defined 
by  the  bracketted  terms  of  equations  (B-50),  (B-51).  and 
(B-52). 


Appendix  C 


Four  Point  Finite  Difference  Approximations  For  Evaluating 

Wall  Derivatives 

The  wall  derivatives,  required  to  evaluate  the  heat 
transfer  rate  and  skin  friction  coefficient,  are  evaluated 
by  use  of  a four  point  finite  difference  scheme.  The  four 
point  finite  difference  approximation  for  the  derivative  of 
a variable  H with  respect  to  the  transformed  coordinate 
is  found  in  a similar  manner  to  that  used  in  the  three  point 
finite  difference  method.  However,  the  four  point  finite 
difference  method  is  more  accurate  since  the  truncated  Taylor 
series,  used  in  the  development,  retain  third  order  terms.  To 
find  the  four  point  finite  difference  expression  for  at 

f station  of  the  grid  in  Fig.  3,  the  following  three 
truncated  Taylor  series  expansions  are  used. 


H 


■*!,  2 


Hi*,’ 


H •*«,  i + i)  H 
3*1 


H«*v  ■+■  ££ 


ATl,  + ! ilH  A*?*  +1 


l 


9 av 


(c-i) 


(fl7?.+A?0  + j.  alE 
2 w 


+ t O. 
9 ao-i5 


-m.i 


(G-2) 


H~,  ‘ 


H/m(  | ■+ 

a* 


+1  (A7*.  +^t  +0*1) 


4_L 

9 3V 


(A*i,  +A>rf+G?i|) 


(C-3) 
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By  introducing  the  definition  of  given  by  equation  (54) 

and  substituting  the  value  of  from  equation  (C-l)  into 

equations  (C-2)  and  (C-3).  the  following  two  expressions  are 
obtained. 


Eliminating 
expression  for 


between  equations  (C-4)  and  (C-5),  the 


dH 

an 


'TH.l 


is  then  found  to  be 


•+I-U,* 


It  K t K' 
K2 


+ H 


3 


V 

(l  +■  K 


~(  l + K tV) 
( l+O  K* 


( C— 6) 


Equation  (C-6)  is  the  four  point  finite  difference  expression 
used  to  evaluate  wall  derivatives. 


Appendix  D 

Computer  Program  Listin, 


ooooonooooooooooooo  oonoooooo 


THIS  COMPUTER  PROGRAM  SOLVES  FOR 
BOUNDARY  LAYER  PARAMETERS  OYER 
EITHER  TWO-DIMENSIONAL  OR 
AXISYMMETRIC  SURFACES 


PROGRAM  TtRACt  (IN  PU  T , OUTPUT , TAPE5=INP'JT  , T APE6=0UTPUT> 

COMMON  G,  p?,  RE  Y , XMINr,  OMEGA,  BO,  TW,  P10,  TIT,  RIO,  VIS10,  TE, 

1 PF,RE,»r , VISINP.SU,  EPS, OS, OYW, SI, ERROR, TC,T A,  I F0GE,IEN01, INTACT, 

2 PRT,XXK,RTRX,XLAMtVARP:*T  ,XINTER,SE!10,ICHS(8)  , I 3RN ( 9 ) * E 0 (200)  , 

3 EN(20n)  ,E°(?00) ,ETO(?on) ,ETN(200) ,ETP(?00) ,F0 (?00) ,FN(20  0) , J2DA, 

4 PP(200),TN(?00),T0(200),XNN(200),YN(200),V0(200),Y3(200),TO(200), 
501(200)  ,T»(  20  0)  , 03(20  0)  , T2JP(200> , T?  JO(  200)  , T2JM  (200 ) 

COMMON/CONST/  R( 70)  , 7C*9)  ,XP  (30),  A(30),  5(30),  0(30),  0(30),  E<  30), 
1F(30),CPC»0)  ,OP(3  0),INIJM 

DIMENSION  Y(200)  ,A1(200, 3) ,A2(200,  3)  ,13  (200,3) ,01(200,3)  , 

1 R2(20n,-»)  ,33(2(1  0,3)  ,01  ( 20  0,3)  ,C2(  20  0, 3)  ,03(20  0,3) 

110  0 FOPMAT ( 1H0,12X,3HS/L,15X,2HCP,15X,6HP/DINP> 

1101  PORMATd  X,3(4X,P1S.9)  ) 

2000  F ORM A T (IX,*PRQFILE  FAILED  TO  RELAX  AT  M = *,I5> 

8002  PORMAT  (5-15,9) 

8003  FORMAT  (1013) 

9002  P OPM AT (1H1,47X*T NTERACTTNG  SOUNDARY  LAYER  SOLUTION*) 

9003  FORMAT  ( 7 HOG AMM A =F3. 3 , 4 4 PR=F6.3,54  MFS=Fo.3,74  R-Y-S=E10.4, 94  T-S 
1 (K)  = P7.  1,11H  HO  = TW/T10  = P6.4, 5H  EPS  = P3.5) 

9004  FORMAT (5H0p10=,E10#4,7H  RHOl 0= , E 10, 4 , 54  T10=, E 1 0, 4, 74  VIS10=,E10.4 
1,4H  SI=,r10»4) 

9005  FORMAT ( 7HOOME" A= , F7,  4,  2X  * 6HPRT  = , c7 • 4, 2X, 7HBT  RX  = ,F7,4) 

9019  F ORMAT ( 1 0X»*WIT4  INTERMI TTENCY  CORRECTION*) 

9020  FORMATdOX,*  WITHOUT  INTE PMITTENCY  CORRECTION*) 

9021  FORMAT  d OX, *T WO -DIMENSIONAL  SOUNDARY  LAYER*)  ! 

9022  FDR«AT(lOX,*AXISYMETRTC',L  SOUND  ARY  LAYER*) 

9096  PORMAT (10X, "WITHOUT  TRANSVERSE  CURVATJRE  EFFECT”) 

9087  PORMATd  OX, "WITH  TRANS  VRSE  CURVATURE  EFFECT”! 

IIIIITTIIIIIITIIIIIIIIIIIIIITIIIIIIIIIIITTIIIimiIITTIIIIItIIIIIIIIIII 


INPUT  INITI AL  CONDITION^ 


G^SPECIPIC  HEAT  RATIO 
PRsSTATTD  PRfiNDTL  NUM3E? 

XMINP=PRrE  STREAM  MAC*  NUMBER 

REY=EYNDLD5  NUMBER  3A5E9  ON  SCALING  LENGTH  AND  FREE  STREAM 
CONDITIONS 
TA=FRFE  STRrftH  STATIC  TEMPERATURE 

OS=INTEGOATION  STEP  SI7-  FOR  STRFAMNISE  COORDINATE 
ST*SURFACE  POINT  AT  WHICH  NUMERICAL  MARCHING  3ROOEO JRE  IS  TO  BEGIN 
OMEGA=0  (TP  VISCOSITY  IS  BASED  ON  SUTHERLANDS  LAW) 

1 ( Tr  VISCOSITY  IS  BASED  ON  POWER  LAW) 


FRROR=AOOEBT  ABLE  ERRO»  TN  INTEGRATION  SCHEME 
XXK=EXPonENT  eOR  NORMAL  COORDINATE  GRID  STRETCHING 
BO=RATIO  OP  SURFACE  TEM3rPATURE  TO  STAGNATION  T EMPE  RAT  JRE 
BTRX=nicT ANCE  ALONG  SURFACE  TO  BEGINNING  OF  TRANSITION  FROM 
LAMINAR  TO  TURBULENT  FLOW 
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PrTsTHRoULENT  ORANDTL  NUMBER 
XTNTEPrTMTE^ITTENCY  FACTOR 

DYW=MIUTMUM  STEP  SIZE  FDR  NORMAL  COORDINATE  IN  TRANSFORMED  PLANE 
TEOGE  = TOTAL  NUMBER  OF  POINTS  IN  NORMAL  COORDINATE  DIRECTION 
TDIFF=0  (FOR  7HREF  POINr  FINITE  DIFFERENCE  SCHEME) 

1 (F03  TWO  POINT  CTMITE  DIFFERENCE  SCHEME) 

I"ND1  = tot^l  N'JMOER  OF  INTEGRATION  STE’S  IN  STREAMWI5E  DIRECTION 
MSP=OEST RED  INTERVAL  IN  NUMBER  CF  STE3S  IN  STREAMWISE  DIRECTION 
J?OA=0  (TWO  DIMENSIONAL  SURFACE) 

1 ( AXT5YMMETpIC  ^O^Y) 

I NUM= NUMBER  0F  PRESSURE  DATA  °OINTS  ENTERED 
IPRrS=0  ( FOR  NO  PPESSJPE  GRADIENT) 

1 (IF  D5ESSURE  GRADIENT  EXISTS) 


ITRAN=3  (FOR  NO  TRANS\/E?SC  CURVATURE  EFFECT) 
i ( FDP  TRANSVERSE  CURVATURE  EJECTS) 

ICON=0  ( TF  BODY  GEOMET»v  IS  TO  BE  COM3  JTED  FROM  PRESSURE  DATA 
L DC A T IONS ) 

1 <IF  BODY  GEOMETRY  IS  COMPUTED  rR OM  SURFACE  EDUATIONS) 
ICHS=STB  TION  NUMBER  WME3E  STREAMWISE  STEP  SIZE  IS  TO  BE  ODJB.EO 
I PRN  = NUMB"R  Dc  STEPS  WHrR E DESIRED  WHERE  DETAILED  PRINTOUTS  DF 
PpnFIL£5  xs  DESIRED 

7=CFNTEPLINE  LOCATION  0F  PRESSURE  OF  PRESSURE  DATA  ENTRIES 
R=RAOTll?  OF  PRESSURE  D A T A ENTRIES 

C D=SURFA CE  3 RES SURE  DIVIDED  BY  FREE  STREAM  PRESSURE 
R,BTRX,AND  Z ARE  scaleo  BY  SAME  LENGTH  AS  USED  IN  REY 
P 1P?=RAT ID  OP  TOTAL  PRESSURE  IN  FRONT  OP  SHOCK  TO  TOTAL  PRESSURE 
BEHIND  SHOCK  (IP  F3EE  STREAM  CONDITIONS  ARE  INPUT  AS  THOSE 
BPHTNO  SHOCK  THEM  P1P?=1  ) 

IITIIITIIIIITITIIIIIIIIITIIITIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII 
PrAO<G,800?>  G,  PR,  XMINP,  REY,  TA 
READ (G ,°  0 0?)  DS,  SI,  OM" G A , ERROR,  XX< 

RffAD (G , 8 0 0?)  BO , BTRX , 0=1,  XINTER,  DYW 

READ (G , 80  OB)  IEDGE,INTACT,IDIFF,IEND1,MSP,J2DA,INJM,IPRES,ITRAN, 

II  CON 

READ (G , « 0 OS)  (ICHS(I),  T = 1,  8) 

RPAD(G,«nOB)  (IPRN(T),  I = 1,  9) 

xlam=.g»btrx 

I F( IPRFR , EQ#  0 ) GO  TO  79 
IF(TNUM»  r9, 0 ) GO  TD  20 
REAO (G ,800?)  DPMAX,P13? 

READ (G,  8 0 0?)  ( CP ( I J)  , I J= 1 , INUM) 

79  CONTINUE 

IF(INUM.ED.O)  GO  TD  20 
IF(ICON.  FD.  1)  GO  TO  4 
READ ( 9,901?)  (R(IJ),IJ  = 1 , INUM) 

4 CONTINUE 

REAP(G,9002)  (Z(IJ) ,TJ=1, INUM) 

WRITE (G, 7 til ) INUM 

7111  F ORMAT (1  OX,  "ALL  ",I?,-  DATA  PTS 
Ic( J2DA»  ED.  0)  GO  TO  81 
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THIS  SUBROUTINE  COMPUTES  RODY  GEOMETRY  FROM  A FINITE  NUMBER  OF 
LOCATIONS 

CALL  CONST  (ICON* I REF) 

*********************************************************************** 
81  CONTINUE 

I F( J?Q  A, SO, 1 ) GO  TO  85 
00  83  KT=  l, INUM 
X°(KT) =7 ( KT) 

83  CONTINUE 
86  CONTINUE 

IPCIPRES.EO. 8)  GO  TO  20 
WRITS(6, 1100) 

XMSO=XMTNp*XMINF 

*********************************************************************** 


THIS  SURPOUTINECOMPUTES  LFAOING  AND  TRAILING  EDGE  PRESSURE 
GRADIENTS 

CALL  SMTRt>R(  BTRX,0PMAX,G,XMSQ) 

*********************************************************************** 

COMPUTE  NONDIMENSIONALI^ING  QUANTITIES 

20  11=  1.  ♦ (G  - 1. ) /2.*XMINF**2 

P10  = (l./(S*XMINF**2)) * ( Z l*  * (G/ (G-l • )) ) 

T 10  = (1 • / ( ( G - l.)*XMI'|F**2))*Zl 
R10  = G*P10/(T10*(G  - 1.  ) ) 

TINF  = i-10/Zl 
TW  = BO*T10 

I F (OMEGA  .ED.  0.)  GO  TO  101 
VTS10  = T 10**OMEGA 

E°S  = ( ( ( G - 1 • ) *XMINF**  2)  ** (OMEGA /2. ) ) /SORT (REY) 

VISINF  = TINF**OMEGA 
GO  TO  10? 

101  TC=198.6/((G-l. )*XMINF»*?*TA) 

VIS10  = (T10**1.5)*(l.  * TC) /(T10+TD) 

E°S  = ((  ( (l.  H198.6/TA) ) * ( ( (G  - 1.  ) * XMI NF*  *2)  * * 1 . 5)  ) / ( ( ( G - t.)*X 
1MINF**?)* (193.  6/TA)))/REY)**.5 
VISINF  = (TINF**1.5)M1.  + TC) / (TINFfTC ) 

102  SU=198. 6 

OUTPUT  INITIAL  CONDITIONS 
HRITE<6,9002) 

WRITE (6,  9003)  G,  PR,  XMINF,  REY,  TA , 3D,  E»S 
WRITE (6,000!*)  P10,  RIO,  T10,  VIS10,  SI 
WRITE (6,9005)  OMEGA , PRT , 9TRX 
IF(XINTrR,ED.l. ) WRITE (R , 9019) 

IFCXINTFP.EO.l.)  WRITE (5,9019) 

Ic(XINTcR.ED.O.  ) WRIT" (3,9020) 

I F( J2QA , CQ.  0 ) WRITE (6, 9^21) 

IF(J2DA.NE.  0)  WRITE ( 6, 90 22) 

IF(ITRAN.EQ.O)  WRITE ( 5 , 3 0 86) 

IF (ITRAM,  "Q.  1)  WRITE ( 5, 3 0 87) 
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C INPUT  TNITI A L PROFILE 

C 

12  M ST ART  s'* 

C INITIAL  TT i THE  STREAHWIGF  LOCATION 

S=SI 

if<j?oa.eo. t>  s=o.o 

0S2=DS1=DS 

nx2ns=nxms=nxos=o. 

sspo=i. 

c initialise  the  streahwisf  location 

Y(i)=o.n 
no  201  LL=2»  200 
OY=XXK*MLL-2) *OYM 
201  Y(LL)=Y(LL-l)+OY 

no  700  LL  = 1,  200 
ni<LL)=n?(L->=OMLL)=*NN<LL)  =0. 

VP(LL) =VM(L.) =VO(LL> =-Y(LL> 

FB(LL)=cO(LL)=pN(LL)=TO(LL)=TN(LL)=TO{LL)=EP(LL)=EO(LL)=EH(LL) 


1 ETP(LL) 

="TD(LL)=FTN(LL)=1.0 

T?JP(LL> 

= T2  JO 

(LL) =T2JN(LL> =1.0 

700 

CONTIN!lP 

no  701  J 

= It 

200 

no  701  T 

= It 

S 

701 

A1(J,T)= 

A ? ( J t 

t)  =A?< J,I)  =R1 (J,I ) = R2(J,I) =OT(J,I>  =01 (J*I) 

i =C2  u,  n =:<  ( J,  I)=0. 

«»REF  = G*XMTNr**2 

TPFF  = (0  - 1. ) *XMINF**2 

RCON=0. 0 

7C0N=n.n 

c 

C INITIALS*:  COUNTERS 

c 

NC=i 

KEN0=TN'JM-2 

IFCICON, rO.  1 ) <ENn=IREF 
ICOUN=MSTART 
I 0=IE0GE 
IG=1 

Io=i 

INOCH=0 
ITCNT1  = 1 
IIN=0 
C 

C ***  OFGIN  FIRST-OROER-TRIOIAGONAL  MATRIX  SOLJTION  *»* 
C 

00  115  m=mstart,ienoi 
IF<J?OA.SO.O>  GO  TO  43 
46  CONTINUE 

IF(NC*  EO«  KEMO)  GO  TO  47 
IFCICON.FO.OI  ZREF=Z<NC*1> 

IFCIC0N.E0.1)  ZpEF=P ( NO) 

IF(ZC0N»GE«7REP)  NC=NC+t 
C 0RO7=SL0P"  OF  TANGENT  TO  SURFACE  Op  OOOY 

C ALPHA=AMOLF  3ETWEFN  TANGENT  TO  SURFACE  AND  CENTERLINE 


81 


oooooo  ao  ooooo  o no 


RCON=RADT'IS  0P  AXISYMMETRIC  BODY 
ZCON=CFNTERLTNE  0ISTANCr 

S=SURFACE  COORDINATE  IN  STREAM  WISH  DIRECTION 

47  OROZ=R(NC)+?»*C  (NC)  * ZOOMIT.*  0(NC)*ZC0N**24-4.*E(  NO  *7C0N**3 
ALPHA=  AT  AN(QROZ) 

Ic (S#LT»  ST)  DZ= (PS2/10. ) * COS (ALPHA) 

IF(S.LT.ST)  ZC0N  = 7C0N«-07 

IF(S.LT,  FT)  RC0N  = A (NC)  + 0 ( NC)  *7C0N+C  (NO  *7C0N**  ? + r)( NO  * ZCON* *3*E  (NO 
1) *ZC0N**4 

IP(S.LT.RT)  s=S+DS2/10. 

IF(S.LT.ST)  30  TO  45 
0T=0S?*C0S(ALPHA) 

ZCON^ZCDM+UZ 

PCON=A(NO«-B(NC)  *ZCON4-C(NC)*ZCON**2fO(MC>  *ZC0N*  *3  + E ( NC ) * ZCON*  * 4 

48  CONTINUE 

I F (H. FO.  ^ST ART)  MP=MSTART 
I F (M, EQ»  T SN01 ) MP=M 
IF(M.FO,  (M/MS°) *MSP)  MP=M 

s =s*os2 

0X2DS  = 0X1  OS 
0X10S  = OXOS 

*************************************  ********************************** 

THIS  SUBROUTINE  COMPUTES  LOCAL  PRESSURE ,PRESSJRE  GRADIENT, 

ANO  TEMPERATURE 

CALL  PR‘SSM(S,XMINP, G,PBS1,0P8G1,TETN-, XME , IRRES , PI  32) 

*****  ****************************************************************** 

PE  = P8G1/PREF 

PP  = 0PBG1/PP.P 

TE  = TETNB/TREF 

UE  = SORT ( 2,  * ( T 1 0 - TE)) 

RE=G*PF/( (G-l.  0)  *TE) 

TR=SU/(TETN-*TA) 

IF  (OMEGA ) 642,676,642 
642  XNUE=TE**0ME3A 
G0T0688 

676  XNUE=TE**1,  5M1.  +198.  5/  (T A *TREF) )/  ( TEU  98.  5/  (T  A *T®E  r) ) 

688  CONTINUE 

COMPUTE  LOCAL  XI  ANO  STEP  LENGTHS 

OXOSsOEPIVATIVE  OF  TRANSFORMED  SURFACE  COORDINATE  WRT 
SURFACE  COORDINATE 

oxos=re*ue*xnue 

IP(J20A.N".0)  0XDS=OXDS<,RCON**2 
I F(M,  FO,  ?)  DX10S=DX2DS=DX0S 

0X2=,5»o«?9*(  (i.  4-QS2/0S1)  *DX1DS+OS1*DXOS/ (OS  1 + 037)  -OS2*OS2*DX?OS/ 

1 (OSl*  (HS14-DS2)  ) ) 

REYNOE=PE«lUE*S/XNUE 

reyext=p-y*visinf*reyno' 

I F(M»  EO,  2)  DX 1 = 0X2 
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X=TRANSFORME0  STREAMWISr  COORDINATE 

IF(M.EQ.?)  X=0X9S*SI 
X=X+0X2 

COMPUTE  STEP  LENGTH  F'JNOTIONS 

Yl=2.« (0X1*2.  *0X2)/ (0X14-0X2) 
IF<niPf  . ED,  1)  Y1  = 2. 

Y 2s  ( (0X1  + 0X2)  /0X1)*2.  0 

Y 3=  (OX 2* 0X2/  <0X1 *(0X1* OX?)  ) ) *2. 0 

Y4=<0X1+9X?) /OX  1 

Y5=0X2/DXt 

TWTF  = TW/TE 

COMPUT"  ALOHA,  OcTA,  ANO  LAM90A 

0UF0X=-pp/(RE*UE*PX9S) 

X AL=UE¥,,c /T E 
XOE=2.  0*X*PUEOX/UE 

6908  L“NGTH=tctig£ 


ASSIGN  THE  MATRIX  ELEMENTS  FOR  THE  FINITE  DIFFERENCE  EOLATIONS 

CALL  ELMA  TX ( M , 0X2, X , X" L , X BE , TR  ,1 DI « , Y 1 , Y2, Y 3, Y4 , Y5, T WTE, I TCNT1 , 
1 ITRAN,A1,A2,  AS,  (31, 92, 37,  Cl,  C2,C3> 


MATRIX  TN VERSION,  SOLVE  FOR  F,  THETEA  AND  V 

CALL  MATr9N3(rP,TP,VP,Dl, 02,D3,A1,91,31,A2,92,C2,A3,93»C3»3,lENGTH 

1 ,200) 

C 

C *****  *********************************  *********  *****************  ******* 

I TCNT 1=TT0NT1+1 
N=IE0GEM 

OY=OYH*XXX*» (TEOGE-2) 

VK=(VD(T-0GE)/(XXKM1.M,/XXK)  ) -V°  ( I EDGE-1)  * ( 1,  -1.  / XXK)  *XXK- 
1 VP(IE0G--2)*XXK/(l.fl./XXt<)  >/0Y 
OY=OYW*XXK** (IEOGE-1) 

K0N=N+5 
0065T0=M»  KOV 
OY=OYW*XXK** (N-l) +DY 
FP(IO) !fP{Ill:1.0 
65  VP(IO)=VP(I- 3GE-1H-V<*0Y 
C INITIATION  OF  SIMILAR  SOLUTIONS 

IFCM.EO,?)  SO  TO  8020 
GO  TO  8018 
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8020  no  8019  T=1 » <0N 

VQ(I)=VM(T)=V°(I) 

F0(I)=rN(T)  =PP(T) 

8019  T0(I)=TM(I)  =T»(t) 

c initiation  simi la»  solutions 

8018  IP=IE0Gr  + l 

l)  ANT  -MCTT*  PROFILES  ITERATIONS 
T AU?=(FP(?)  -e-p(l)  ) / 0 Y W 
I F(ITOMT 1 , E0«  2)  TAU1=10. *TAU2 
Rt1?  = TAJ1/TA'J?-1. 

T AUl  = T A'JR 

I r (TTCMT 1 .IE.  100)  GO  TO  7005 
WRITE (8»  2 003 ) M 
CALL  EXrr 

7005  IF(A9SC»Ti?)  .GT,ER=»OR)  GO  TO  699e 
U AND  THETA  o^O^ILES  ITrRATIONS 


COMPUTE  9LT,  nOT (CELT A c T AR)  AN0  BMT ( THETA) 

55  C 0=TD ( 1) 

TPI=0. 

RLT  = 9L  OT  = 9L  N T = 0 • 

XNN(1)=0. 

OO  57  N=7,K9N 
OY=OYW*XX<** (V-Z) 

CX=TP  ( N> 

TPI=TPH-.5*0Y*  (CO+CX) 

cn=cx 

XNN=COOPOINATE  NORMAL  TO  SURFACE 

XNN (N) =r®S*  TDT*  SOpT  ( 2,  *X)  /<RE*UE) 

IF(JPOA.N-.O)  XNN  (N)  =XNN  ( N)  / RCON 

IF(J2DA.E0.  1)  T?JP(N) =1.  +2.*S0RT(2.*X>  *EpS*C0S ( ALPHA ) *TP I/( RE*UE 
1*RC0N**2> 

IF(TTRAM,  NE.  0)  XNN(N)  = ( ® CON/COS  (AL ®H A > ) * (SORT ( T 2 J®< M)  ) -1 . ) 

9L0T  = 9LOT  + (2,-PP  (N)  /Tp(N)-FPCN-l)  /TP(N-i)  ) *CXNN(N)  -XNM(N-l)  ) /2. 
9LHT  = PLMT  + (-P(N)  * (1. -rP  ( N) ) /TP  ( N)«-F°  ( N-  1 ) * Cl. -r  P(N- 1 ) ) / T P (N-l ) ) 

1 *(XNN(N)-XNN(N-l))/2. 

TF(BLT.GT.O.)  GO  TO  57 

IF(FP(N)  .GE.  0.995)  9LT  = XNN  (N)  - (FP(  N)  995)  * (XNN  (N>  -XNNCN-l)  ) 

1 /(FP(N)-FP(N-l)  ) 

57  CONTINUE 

COMPUTE  9LT,  90T ( OELT A STAR)  ANO  BMT (THETA) 

COMPUTE  THE  E09 Y VISCDSTTY  COEFFICIENT 
IF(S.LF.BTRX)  GO  TO  55 


THIS  SUBROUTINE  COMPITES  EDOY  VISCOSITIES 

CALL  REYSTR  <KONt TR, X , TREF, XNUE , X8E, S, I TCNT1 , 
1RCON, AL°HA,  ITRAN) 
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58  ITCNT 1 = 1 

ASSESMFNT  OF  GRID  PONI T'  IN  ETA 

IF(INOCM)  71,  71,  772 
71  CONTINUE 

I p (M  - 70)  732,  732,  ?2 

7?  IF  (AOS  (ro  (I"0^“-15)-FD(rrnGE-16)  )-0,  0 001)  73,77,7!, 

73  IC(ABS(TO(IEOGE-1E)  - T°  CIEOGE-16)  > - ,0001)  732,  732,  74 

74  IFOGF=T-nr,-f  i 
TQ=IQ*1 

OY=OYW*XXK**  <TEPGE-2) 

Y(IEOGF)  = Y (TrOGE-l)  ♦ OY 
732  IQ  s IQ  - 1 

ASSESMENT  OF  GRID  PONITS  IN  ETA 


THIS  SU°RQUT I ME  COMPUTE'  HALL  STRESS  AND  MEAT  TRANSFER  RATES 

CALL  CPSTNO  ( T P,  XNUE , X,  S , X9E,  M , BL  DT  , BL  MT  , BLT  , ®3Gl , 0P8G1  , REYEXT  , 

1 XME»MP,PCON,ITRAN,7CON) 

*******#*****************»***¥*****************»**********«,******»♦***» 


shift  profiles  back  one  xi  station 

NN  = IQ  ♦ 5 
OO  118  N=1,NN 
FN(N)=FQ(N) 

FO  (N)  = co  ( N) 

TN(N)=TO(M> 

TQ(N)  =TP(N) 

VN(N) = VO (N) 

VO(N) = Vp (N) 

ETNC  N) =ETQ( N) 

ETO(N)=ETp<N) 

FN(N)=EO(N) 

EO  CM)  = EP  ( N) 

T ?JN (N)=T,JO(N) 

T2JO(N)=T2Jp(N) 

118  CONTINUr 
OXl=OX2 
0S1=0S  2 

Ie(M+l-rCHS(IG) ) 114,117,114 

113  OS2=2.0*QS1 
IG  = TG*1 
INOCH  = 1 

IF  CM.  S3.IEN01)  GO  TO  27 7 
GO  T0  111 

114  OS?=OSt 
INOCH  = 0 


oo  o a n n o o o n n a o o n n 


IF  <M,  E9. IEM01)  GO  TO  2*7 
GO  TO  ill 
237  TIN  = 1 

************************#*»******♦•*********»**♦********¥***»********** 

THIS  SUB°OUTINE  ALONG  WITH  SUBROUTINE  CFSTNO  IS  USED  TO  OUTPJT 
ALL  OATS 

111  CALL  PRMOHS  ( ICOUN , I P , IG , 1 0 , MSTA RT , IIN, M, S , f ,«L T ) 


115  CONTTNIF 
rNO 


///////////////////////////////////////// nt /tr  ft  nt  tt/t/t/tt  /ft  ft  rtfft 
SUBPOUTIN'  P*ESSM(S,XM,G,»,OPOX,T,YH,IPRES,P13?> 
///////////////////////////////// //  mini  a //////////////////// /////// 

COMMON/*  OUST  / R(30)  ,7<3r>)  ,XP(30>  ,A(30>,  O(30),C<39>  ,0*30)  ,r<30)  , 

IF  (30)  ,0^(30)  ,OD  (30)  ,INU‘4 

100  FORMAT(  F*,*WARNING....,'ALCULATION  IS  OUTSIOE  OF  THE  PRESCRI3E0  PR 
1F<5SURF  OATA,  S IS  LESS  THAN  XP<1)*> 

200  F0RMAT<  5X,*WARNING....~ALCULATION  IS  OUTSIOE  0r  THE  PRESCRI3EO  PR 
1ESSURE  TATA,  S TS  GREAT'S  THAN  XP(ENO)*) 

300  FORMAT (1X,5E15. 9) 

IR=n 

IPM1=INUM-1 

IFdPOFR.RO,  o)  GO  TO  <*0 
00  20  1=1 » IMUM 
IF(S.LT. XP(1) ) WRITE { 3 » 1 0 0) 

IF(S»GT,X°(IN,JM)  ) WRITE  ( 6 » 20  0 ) 

IF (S«LF » Xn( 1 ) ) IR  = 1 
IF(IR.NF.O)  GO  TO  30 
IF(S.GE.X°(I»M1)  ) IR=  I N'J M 
IF(IR.NF.O)  GO  TO  30 

IF((S. GF.XP(I)) .AMO. (S.LT.XP (I+l)) ) IR=I 
IF(IR.Fn.n)  GO  TO  20 
C SEEKING  THE  BEST  FIT 

Rs=<s-xp(i)>/<xp(m>-xp<i>> 

IF (PS»  GT»  0»  5 ) 10=1*1 
C SEEKING  THE  BEST  FIT 

I F ( IR , Nr%  0)  GO  TO  30 
20  CONTINIJF 

30  IFCTR.GT.IPMl)  IR=IPM1 

TRP=IR  >1 
IRM=IR-1 

IFCIR.FO.  1)  IPM  = IR4-2 

C COMPUTE  THE  CUBIC  S°LTN"  COEFFICENTS 

Xl=(X°(TRP)fXP(IRM)-2.0*XP(IR) ) ♦ (XP( IRP) -XP( IRM) ) 
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X?=(XP(TRM)-XP(IR)>*(X0(TRM>-XP(IR)  ) 

X3=(XP(TPP)-XP(IR))*  (XP(IRP)-XP(IR) ) 

X4=XO( IPM)-XP(IRP) 

X*>=XP(I°M)-XD(  I R ) 

X6=XP(T°P)-X°(IR) 

nETS=X5*X5*XL 

02=(0P(T3)*XH-DP(IRD)  *X2-OP(IRM)*X3)/OETS 
C3=(0P(TR) *X4-0P(IRP) *X?*0P(IRM)*X5) /OETS 
C COMPUTE  THE  CU3IC  SPLINE  COEFFICENTS 

nxp=s-xo(TP) 

oxp2=oxp**2 

OXPFsOXP/pO. 

TPDX1=0B(IR) 

P=CP  ( I °) 

no  10  1=1,20 

x=l*oxpr 

x?=x*x 

0°0X2=  O'*  ( TR)  fC2*X+C3*X2 
p=p+o.p*  (opdxh>opox?)*oxpf 
10  o°oxi=o»nx2 

opox=op(I?)  ♦ C2*nyp«-co»oxP2 

T=(P*P1P2)** (( 0-1.0 )/G) 

Y M=S0RT (2.  0* ( (2, 0-MS-1. 0)  *XM*XM) /(  2.  0 * T ) -1. 0)  / ( G-l . 3 ) ) 

WRITE(G,300)  S,P,DPOX,T,  YM 
GO  TO  F0 
40  P=l.  0 

DPOX=0. 

T=  CP»P1P2) ** ( (G-1,0) /G> 

YM=SORT(?.0* ( (2, 0+( G-1,0) *XM*XM)/(2» 0 * T ) -1. 0) / ( G-l. 0 > > 

50  RETURN 
END 

/////////////////////////////////////////////////////'///////////////// 
SUBROUTINE  5MTHPR (BTRX , OPMAX ,G , XMS0) 

// / // / /// // /// /// ///  / r n / // / //  / // // / / // // ///  nr  urn  t rr  //  nr  r /ft  n //  rt  / 

COMMON/? ONST/  R ( 30) , Z ( 30 ) , XP (3 0) , A (3 0 ) , B ( 30) , Z ( TO ) , D (TO) , E ( 30 ) , 

IP  (30) ,CP(T0) , OP ( 30) , INUM 

100  FORMAT (1 X, 'FIRST  OP  OATS  POINT  YIELDS  ADVERSE  PRESSJRE  GRADIENT  TO 
10  STEEP  POR  CALCULATION  TO  CONTINUE*) 

200  FORMAT (1H0 , 1 1 X , 3HS/L , 15X , 2HCP, 11X, 6HP/» INr, 14X , 4HOP0X) 

300  FORMAT (IX, 4(VX,E15.R) ) 

DPTOL  = OPMAX*l,  01 

C COMPUTE  THE  TRAILING  EDGE  DPOX 

TPM1=TN'JM-1 
I °M2=I MUM-2 
OXi=XP(TOMl)  -XP(IN’JM) 

OX2=XP(TPM2)  -XP  (INUM) 

0X12=0X1*0X1 

0X22=0X2*0X2 

OP(INUM)  = (OP(IPM2)*OX12-CP(IPM1)*OX?2-:P(INUM) * ( 0X1 2-DX22)  ) / 

1 (OX1*OX9*(OX1-OX2)) 

C COMPUTE  THE  TRAILING  EDGE  OPOX 
10  I MAX=  0 
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C COMPUTE  THE  LEADING  EDGE  OPOX 

oxi=xp(?)-x»(1) 

nx?=XP(T»-x»(i) 

0X12=0X1 *0X1 
nX2?=0X'**DX2 

OP (1 ) = (P° (X ) *0X1 ?-C° (2)*DX22-CP(1) * (0X12-0X22) ) / (OX  1*0X2* (0 XI- 
1 0X2)) 

10(00(1) . OT.OOMflX)  WRIT- (6, 100) 

TP(0O(l).  r»T.  OOMAX)  CALL  EXIT 
C COMPUTE  the  LEADING  EDGE  OPOX 

no  20  1=2 , T 3v11 
I *11=1- 1 

tpi=tm 

OX1=XP(TM1)-XO(I) 
nX2=X»(TDl) -X° ( T ) 

0X12=0X1 *0X1 
0X22=0X0*0X2 

OP(I)=  (r°  (I PI)  *0X12-CD(TM1)*0X22-CP(I)*  (0X12-DX2*) ) M0X1  *0X2* 

1  (OXl-nx?)) 

20  IF((OP(T).r,T,OPTOL).ANO.  (XP(I)  .LE.9TRX)  ) IMAX=I 
I p (I MA X»  r 0.  0 ) GO  TO  50 

C SMOOTH T^G  THE  CP  OATA  IN  THE  LEADING  EDGE  REGION 

I MM1  = 1 MS  X-l 
I MP1  = I Mfi  X*1 
0X1=XP (TMM1 ) -X° (TMAX) 

0X2=X° (TM^l) -Xp ( IMAX) 

0X12=0X1*0X1 

0X22=0X2*0X2 

CP(IMM1)  = (CP (IMP1)*0X12-CP( I MAX) *(0X12-0X22) -OX  1*0X2* (OX  1-0X2) 

1 *0PMAX)/0X22 
GO  TO  10 

C SMOOTHING  THE  Cp  OATA  IN  THE  LEAOING  EDGE  REGION 

50  WRITE  ( ft, 2 00) 

00  30  I=1»INUM 

PC=2.  0*(C°(T)  -1.  0)  / (S*XMSO) 

30  WRITE  (6,100)  XP(I),PC,Cp(T)  ,OP(I) 

RETURN 

END 

/////////////////////////////////// 

SUPROUTTNE  CpSTN0  ( TR, XNUE, X ,S , XBE, i , 9L OT , 9LMT , BLT , PBG1 , DPBG1 , 

1 REYEXT,XME,MP,RCON,ITRAN,ZCON) 

tutrrt/ttrtitmttHtnntiniiiinittuittfutntitwttttnnntint 

COMMON  % PR,  REY,  XMINP,  OMEGA,  3D,  TW,  P10,  TIT,  RIO,  VIS10,  TE, 

1 PE, RE, UE, VIS  INF, SU, EPS, OS,OYW,SI,ERRDR,TC,TA,t  EDGE , IE NOi, INTACT, 

2 ORT,XXX,3TRX,XLAM,VAPPPT,XlNTER,SE®0,ICHS(8),IPRN(0),E0 (200) , 

3 EN(200) ,"p(200) ,ETO(20n) , ETN ( 20 0) , ETP( 20 0) , PD ( 210 ) ,FN(?00> ,J20A, 
l*  FP(20  0)  ,TN(200)  , TO  (200)  , XNN  (2 00  > , VN  ( 20  0)  , VO  (20  0)  , V=>  (20  0 ) , TP ( 20  0)  , 
5 01(20  0)  ,02(200) ,03(200) , T2J°( 20  0) , T2 JD (20  0) , T2 JN(20  0) 

2999  FORMAT (1H1, 10X, 5HS/L  = , “ 15. 8 , 1 OX ,6H  Z/L  =, E15. 8 , 10X, 5HR/L  =,E15.8) 

2000  FORMAT (!H0,1PX,5HS/L  =,E15.8) 

2001  FORMAT (’X,THXME  =, El  3.  8 , 2X , 7HPE  =, E15,8,2X,TH0PPIXF  = ,F15. 8 , 

1 ?X,7HXPE  =,E15.8,2X,7HTW/TE  =,E15.3) 


CIO  oooooo 


"2002  FORMAT <?X,7HRLT  = ,F1E.  3 ,2X  ,7mr\_MT  = ,E15.  8,2K  ,TH«3lOT  =,£15.8, 

1 ?X,7H«?'Y1T  = , E15.8,2X,  7H?5YDT  =,£15.8) 

2003  F0RMAT(M,  7NCPN0  =,£15.  8,2X,7HCF£N0  = , F15.  8 , 2<  , ’HSTNO  =,E15.8, 

1 ?X,*HGTEN0  =,E15.8,2X, 7MREYEXT=,E15. 3) 

TMTE=TP(1 ) 

IFCOMEGA.EO.O, ) GO  TO  3=5 
IF (OMEGA  .£1.  1.)  GO  TO  3551 
XLM1  = 1,/(TMTE**(1.  - 0*EGA)) 

GO  TO  856 

855  XLM1  = <1.0  + TR)  *SORT < TWT r ) / ( TWTE+TR) 

GO  TO  85  5 

8551  XLM1  = 1. 

856  CONTINUE 

Yll=  ( ( ?.  + XXO  * ( 1.  +XXK  + XY»<**2)  U.+XXO  /(  (1.  + XXO  Ml,  +XX<+XXK**2)  ) 
Y12=(1.+XXK+XXK**2> /XXKR*2 
Y13=(1.+XXK+XX'<**?)/(XX<**3*  (1, +XXK) ) 

Y14=l.  /(XX<**3‘  (l.+XX<+YXK**2)  ) 

TAU=XLM1*RE+XNUF*UF*UE*  ( - Y11*FP  ( l)  + Y1 2*  FP ( 2)  -Y13*F»<  3)  +Y14*FP<4)> 
1 /(OYWMORT  (2.  *X)  ) 

OF  = XLM1*RE*XNUE*UE*TEM  Y11*Tp<1)  -Y 12* TP (2)  +Y17*Tp<3)  -Y14*TM4> ) 
1 /(OYWMORT  (?.  *y)  *PR) 

IF(J20A.  NE.  0)  TAU=TAU+R~ON 
IF(J20a.N-,  0)  QS=QS*R:0N 
5TN0  = n. 

IMRO  J = , 1,1  5TN0  = £3S*05/((1.  - OOIMTE  + .5 *U£**2) ) 

STENO  = FTN0/(RE*UE) 

CFNO  = ?• *£3S*T  au 
CFENO  = CrNO/ (RE*UE*UE) 

REYOT=RE YEXT*RLOT 
RrYMT=RrY“XT*OLMT 

selection  OF  THE  OUTPJT 
IF(M.NF.mo)  jo  TO  1000 
SELECTION  0=  THE  OUTPJT 


OUTPUT  STATION  DATA 


IFCJ2OA.E0.1)  GO  TO  457 
WRITE  (6,2100)  S 
457  CONTINUE 

IF(J20A»  rO, 0 ) GO  TO  456 
WRITE ( 6, 2099)  S,ZC0N,PC0N 
456  CONTINUE 

WRITE  (6,2001)  XME,POSl,">POGl,XRE,TWTE 
WRITE (6,  2002)  BLT^LMT^LOT.REYMTjREYOT 
WRITE (6, 200 3)  CFNO,CF£NO, STNO, STENO, REY EXT 
10 00  RETURN 
ENO 

// // / / / / / / / // / / / / t /// / / / i / / // / / / / / / / / ft / / ff f / / r f / ( f t / r r r / r / / / / / / f r r / rr / 

SUPPOUTTNE  ELMATX  ( M,f)X2»X,XAL»X0E»TR,I0IcF»Yl,Y2,Y3»Y4,Y5,TWTE, 

1 ITCNT l,ITRAN,At,A2,A3,'M  ,P2,R3,C1,C2,C3) 

C ////// ////////////////7///////////////////////7Z/////////////////7///// 


89 


occ  OUOC(J 


COMMON  r-,  PR,  RPY,  XMINC,  OMEGA,  93,  T W,  °10,  TIT,  RIO,  VIS10,  TE, 

1 PEf‘?Et0"t'/I?TNP'*SU,£:,S,r>^,0YWfSI,ER,?3?fTCfTA,IF0GE,  IENOl , I NT  ACT, 

2 °RT  ,XXY,9TPX,XLA'SVAPP::>T,XINTER,SEP0,irHS(8)  , I P?N  (9), £0(200)  , 

3 CN(200)  ,EP(200)  ,ETO(20C)  ,lTN(2  00>  ,ET»(?00)  ,F3(210)  ,FN<200)  , J?0A, 

u pp(?on>  ,tmc?oo)  ,To(?on)  ,xnn(?oo),vm(?ooj  , vd(?oo>,vp(200),tp(2oo>, 

‘»01(20  0>,0'>(?0n),03(?10),  T?JP(200),T?JD(200)  , T2JNC203) 

OTM£NSTON  Al(?00,3)  , A?(',00*3)  , A 3 (20  3 , 3)  ,91(200, 3)  ,92  (200,3)  , 

1 03(200, 3)  ,31(200,3) ,D2(?00,3) ,C3(200»3) 

THF  TNNrR  EDGE  BOUNDARY  CONDITION 

HO  8011  1 = 1,3 

8011  Al(l,T)=A’(l,I)=A3Cl,I)=31(i,I>=82<i,I) =93(1,11  =31(1,1)  =C3(1,I) 

1 *C3(1 ,T) =0. 

A 1 (1, 1 ) = 1. 0 
B?(i,n*i,o 

D1(1>  = 0. 

Q?(l) =TWT- 

IFCSEPO#  rD.  0#  > GO  TO  001? 

C3(l,l)=l,0 
03(1)  =0, 

GO  TO  0013 

8012  XL=OX?/(?.0*DYW> 

A3(l,l)=nx2*X*Yl 

C3(l,i)=-2.  *XLM2.+XX<) Ml.  + XXK) 

C3(1,?)=?.*XL*( 1. ♦XXX) /*XK 
C3(1,3)=-?.*XL/ (XX<M1.*XXX) ) 

03(1)  = 0. 

THE  INN"?  EDGE  «OUNOARY  CONDITION 
THE  FIELD  POINTS  EVALUATION 

8013  NMl=I~nGr-i 

DO  8014  *i-2f  NM1 
DY=XXK** (N-l) *OYW 
0YM1=DY/XXK 
XL=OX2/(P.O*OY) 

Y6=?./(t.*0YMl/0Y) 

Y^sDY/nYMl 

Y 8=?,  /((DYMl/DY)*(l.OYMl/0Y)) 

Y 0=2. 7(1.  ♦0Y/DYM1) 

Y 10=1* -DY/OYMl 
SEP=1. 0 

IF(FO(N).LE.O.)  SEP=0, 

IFCITCNT1  . 3 T . 1)  GO  TO  7000 
IP(IOTFr  a ED « 1)  GO  TD  7501 
FMl=Y4*rD (N)  -Y5*FM(N) 

TM1=Y4*T0(N) -Y5*TN(N> 

VM1=Y4*V0(N)-Y5*VN(N) 

IF(SEPO.  ED.O.)  VMl=VO(N) 

EMI  = (Y4*  (EO(N-l)+EO(N)  +EO  ( N+ 1 ) >-Y5*  (E  N (N-l)  f-N(  N)  *EN(N+1>  ) ) 73. 
ETM1  = (Y4MET0(N-1)  «-ETO  ( N)  ♦ ETO  (N+ 1)  ) - Y 5*  ( ETN(  M)l  ETN(  N H-ET  X ( N*1 
1 >)>/3. 


5.'i!T7"T 


T?J=Y4*T?  JOCM>-Y*5*T2JM<M> 

o o to  ’’ooi 

7? 01  Fmi  s F0(N) 

T Ml  = TOCO 

V Ml  = VO(O) 

*=■^1=  (~0(  M- 1>  *20  (N)4-"0(NM))  / 3. 

ETM1  = (ETO  (N-l)  +ET0(N)  «-ErO  (N  + l)  ) /3. 

T 2J=T? J® (M) 

GO  TO  7001 

7000  FM1  = r<MN> 

T Ml  = TP<M> 

V Ml  = V°(M) 

EM1=(F°(M-1)  t-P(N)f£P(NH)  )/3. 

ETM1  = ( FT ° (M- 1)  + ETP<N)  f£-P<N+l>  ) /3. 

t?j=t?j»(m) 

7001  IF (OMEGA  .-O.  0.)  GO  TO  634 
Ic(OMF.G5  .E3.  1.)  GO  TO  6341 
XLM1=1./ (TM1**( 1. -OMEGA) ) 

XLPMt=  (OMEGA-1. )*XLM1/TM1 
G0T0625 

6841  XLM1=1. 

XLPM1=0. 

G0T0625 

684  XLM1=  ( (1.  ♦■TP)  *S0PT(TM1)  / (TM1+TR)  ) 

XL°M1=XLM1*  (TR-TM1)/  ( ?.  * T M1M  TM1+TR)  ) 

625  I F (IT0NT 1 « GT  # 1 ) GO  TO  626 
FY=(Yq*«*0(M*l)  / ?.-riO»FO(N)-Y3*FO(N-l)/2.)/OY 
TYr (YO*TO(M+l) /?.-Y10*TO (N) -Y8*TO(N-l) /2. ) /OY 
EYMl  = (Y9*20(N«-i)  /?.  - Y 1 0 * EMI- Y8  *EO(M-  1 ) /2.) /OY 
ETYM1= (Y9*2T0(N*l)/2. -Y1 0 *ETM1-Y9* ET 0 ( Y-l ) /2. ) /OY 
0T2J=  (YP*T?J3(N|  + d/  2.-Y1  0*T2JO(N)-Y3*T? JO(N-l) /2.) /OY 
GO  TO  6 77 

626  FY=(Y9*rP(Mf  1)  /?.-Y10*cP(M)-Y3*FP(N-l)  /2.)  / DY 
TY=(YO*tP(M  + 1>  /?.-Y10»Tp  (M)-Y8*T°(N-1)  /2.)  /OY 
E7M1= (Y3*"P(N+l>/2.-Y10*EMl-Y8*EP(N-l) /?.) /OY 
ETYM1= (Y9*ETP(N+l)/2.  -Y1 0 *ET M1-Y8*ETP( N-l) /2.  ) / nv 
OT2J=(Yn*-?jD(N+l)/2.-Y10*T2J® (N)-Y3*T2 JP(N-l) /2.) /OY 

627  IF(IOIFF. ^0. l)  GO  TO  7502 
FM2=Y2*F0(M) -Y3*PN(N) 

TM2=Y2*T0 (N) -Y3*TN(N) 

GO  TO  7?05 
7502  F M2  =2.*fo(M) 

T M2  =2.*T0(M) 

7505  CONTINUE 

IF(ITRAM,E3.  0)  T2J=1.0 
IF(ITRAM,EO. 0)  OT2J=0. 

A1(N*1)=Y3*XL*( 2.*XLM1*-mi*T2J/0Y-XLM1*FM1*0T?J-XLM1*T?J 
1*EYMH-Vm1-EM1*T?J*XLPM1*TY) 

Al(N,2)  = -4.*Y^*XL*XLMl*EMi*T2J/0Y-2.  *Y10*XL»  CXLM1*EM1* 
10T2J+XLMt*T?j*EYMl-VMlf  7‘»1*T2J*XLPM1»TY)  -?,*DX2*XOE* 

?FMl*SEP-X*(2,  *Y1*FMI--M?) *s£p 

Al(Nf 3)=XL*(2.*Y6*XLMl*"Ml*T2J/0Y+Y9*(XLMl*EMl*0T2JfXLMl 
1*T2J*EYM1-VM1*EM1*T?J*XLPM1*TY) ) 
B1(N,1)=-Y3*XL*EM1*T2J*XLPM1*FY 
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B1  (N,  ?)=-'*•  *Y1  0*XL*EMl*r?j*xiPMi*PY4-OX?*XBE 
0 t (N  , *)  = Y0*XL*EM1*T2J*XI.  DM1*FY 
C1(N,1)=C1(N,3) =0. 

Ct(N,2)=-DX2*FY 

A?<N,1>=-’.  *Ya*FY»XL*XLHl*ENl*T2 J*XAL 

Ab(N,2>=-4.*Y10*FY*X|_*XIM1*EM1*T?J*X  AL-X*(  Y1*TH1-th?)  *SEP 
A?(N,  3)=*>.*Y9*FY*Xt*X_Ml*EMl*T2J*XAl. 

82(N,1)=(XI*Y3/PR)*(2.*XIM1*ETM1*T2J/0Y-XLM1*ETM1*DT2J- 
lX1_Mi*ETVHl*T2.H-PR*VNl-2.  *TY*XLPM1*ETM1*T?J) 

8?(N,  2)=-4.  *Y7*XL*XL^1*"TM1*T?  J/  (OY*  PR)  -2.  *Y10*  (XL/bR>  * 

1 ( XLMl*"TMl*!3T?J'fXLk11*"TY‘*l*T2J-PR*\Ml>?«*XLPMl*rT'il* 
2T?J*TY)-X*Y1*CH1*SEP 

D2(N,3)=(XL/PR>*  (2.*Y6*XLM1*ETH1*T2J/:)Y+Y3*  (XLH1*ET11* 
10T2vHXLm1*ETYH1*T2J-bR*'/m1  + 2.*TY*XL°M1*ETmi*T2J>  ) 
C'»(N,2)=-DX2*TY 

=0. 

A3(N,i)=A3(N,3>  =0. 

A3(N,2)=0X2«-X*Y1 
83<N,t>=83(N,2> =B3(N,3>=0. 

C3(N,  1)  = -XL*YF 
C3(N,2)=-B.  *XL*Y10 
C3(M,3)=Xl*Y9 

Ot  (N>  = nv?*PY*  (EM1*T?J*XLDM1*TY-VM1)  -P'11**2M0X?»X3Ef  X* 

1 Y 1)  *SFD 

O2(N>=0X?*TY*<XLPMl*ETHi  *T2 J*TY/PR-(/Ml>  ♦0X2*XLN1*EH1*T? J 
l*XAL*FY**'>-X*Yt  *TM1*FM1*  SCP 
03(N)=X*FM2 
8014  CONTINU- 

THE  FIELT  POINTS  EVALUATION 


THE  OUT-R  EDGE  BOUNDARY  CONDITION 
00  8015  1*1,3 

8015  Al(TEOGE,T)=A?(TEnGE,T)=A3(IEnGE,I)=31(IEDGE,I> =32 (I EDGE, I) =3  3 ( 
1 IEDGE ,T) =C1 (IEDGE, I) = 32 ( IE0GE , I )=C3 ( I l 0GE, I) =0. 

Al(IE05i,7)  =i.  o 

B?(IE0GE , 3) =1, 0 

OKTEOGD  *1.  0 

02  (IEDGE)  =1.0 

I F(SE°0»  ED»  0 • ) GO  TO  3016 

XL=DX2/(B.*0YW*XXK*MIEOGE-l)> 

FM2=Y?*-0 (IEDGE >-Y3*FN(TEnGE) 

I F ( IOI Fr , 20, 1 ) FH2=2. *FD(IEDGE) 

A‘?(IEOOE,3)=0X2  + X*Y1 
C 3 (IEDGE,  1)=2.*XXK**3*XL/  (1#  ♦XXK) 

C 3 (IEDGE,  2)  = -2.  *XXK* (1.  «-XXK>  *Xl 
C3(IE0G",3)  = 2.*XXK*XLM?.*XXK+1.)/(1.*XXK> 

03(IEDGE)=X*FH2 
GO  TO  80  1 * 

8016  VH1=V0 (TEDGE) 

IF(ITCNrl»GT • 1)  VM1=VP(TEDGE) 

C3(IEOGr,3)=1.0 
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OOO  DO  DO  DO 


03(IEDGr) =VM1 
8017  rONTINUr 


THE  0MTr  R E3G"  POUNOA’Y  CONDITION 

RETURN 

END 

///////////////////////////////////////// ///  / / / 1 f / r t r / r / 1 1/ ( / f / / /////// 


SUBROUTINE  P°NCHS(ICO'JN,Ip,IG,  IQ ,MST  APT  , IIN  , M , 5 , Y,  BLT) 


///////////////////////////////////////////////////////////// ////////// 
COMMON  G , 3R,  REY,  XMIN",  OMEGA,  90,  rN,  RIO,  Ml,  RIO,  VIS10,  TE, 


1 °E,R",HE,YISINP,S*J,E:»G,pE,OYW,?r,ERROP,TC,TA,IEOGE,  IEN01 , INTACT, 

2 °RT  ,XXY,9T?Y,YLAM,VA=,3',T,XTNTER»SE°0,ICHS(3) , I °RM  ( 9 ) , 10  (20  0 > , 

3 EN(200) ,EP( 200)  ,ETn(?0fl)  ,ETN(  20  0>  ,£T*(210)  , FO ( 200 >, fn ( 20 0)  ,J20A, 

4 cp(?nij)  ,TN(?nO  ) ,T0(200)  , XNN  (200), VN  (20  0)  , VO (200)  , V=>(200>  ,TP(200) 
5Dl(2  00),D’(200),  03(20P>,'r2JP(200),T7J:)(?0  0),T2JN{?00> 


DIMENSION 

Y (200)  ,7(7,16) 

25 

FORMAT  f 1H0, 45X , 

23HPR0FJL r FOR  STATION  S =F14,3) 

40 

FORMAT (°MTN= 

15F3.4  ) 

41 

FORMAT (»H 

ET  A = 

15F0.4  ) 

42 

FORMAT  («M 

Fl  = 

15F3.4  ) 

43 

FORMAT (OH 

T 1 = 

15FS.4  ) 

44 

FORMAT  ( 0 H 

V 1 = 

15F3.2  ) 

46 

FORMAT  (6H 

E0  = 

15p0. 2) 

507 

FORMAT  (PH 

Y/OLT= 

15F3.4  ) 

509 

FORMAT  («M 

RD /ROE 

=15^3.4  ) 

510 

rOPMAT  («H 

MACM  = 

15F3. 4) 

511 

FORMAT ( « M 

DT/POP 

=15F3,4  ) 

512 

FORMAT  (PH 

PT/PE= 

15FS.4  ) 

513 

FORMAT (PH 

H / H E = 

15F3.4  ) 

IFdCOUN- 

TPPN(IP) ) 51, 3?, 51 

* 


OUTPUT  °ROFI LE  DATA 


38  KONT  = 10-1 
J 7=0 

WRITE ( 6,23)  5 
0050J1=1,KONT,1E 
J?=J2+1 
KON=  J?*1  5 

WRITE  (F,40)  (XHN(N) ,N  = J1 ,KON) 
WRITE  (6,41)  (Y <N),N= Jl, YON) 
WRITE  (6,4?)  (PO(N),  N=  J 1 , KON) 

WRITE  (6,43)  (T  O ( N) , M=J1,K0N) 
WRITE  (6,44)  ( VO  ( N)  , N=J1,K0N) 

WRITE(6,46)  (EO(N),N=Jl,KON) 
I=J1-1 

IF(M,EO,  MSTART)  GO  TO  50 
00530JX=1,15 
1*1*1 

2 (1, JX)=XNN(I) /HLT 
Z(2,JX)=F0(I) 

C 7.  (3,  JX  ) = TO ( I ) 
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7(3,JX)=l./T0(I) 

ptppo=  (G-i.o)  *te:*to(I) 

T F (PT  Prn ) 777,777,778 

777  PTPE0=1. 

778  7 (4, JX) =Ur*T (’, JX)7(PTo-n)»».5 

PTPF  = 7 (t»,  JX)  ** ( 4, JX) 

TF(Z(4,JX)-1.  O)  504,504,^05 

504  P-rPc!=  (1,  n«-(((3-1.0)/R,0l  •PTPE)  ) ♦*(  37(3-1.0) ) 

GOTO50F 

505  PT P=>(  ((CM.  0)*PTPE/7. 0)  **  (G/CG-l.  0)  ) ) M C ( 3*  1 . 0 ) / ( ( 2 , 0 * G*PTP£> 
11,0))  ) *M1.0/(G-1.0>  )) 


506 

530 


50 


7 (5, JX)=3Tpr 


(6, JX ) = ptdT'PF/ P10 

(7,JX)s(TE*T0(I)  /(UE*U'_)  + .5*F0(I)*r0(I>)/  (TE/ (UE*'JE>  ♦.  5) 
)NTIHU“ 

?ITH(F,^07)  (7(1, N)  ,N=1,15) 

RITE(6,500)  (7(3, N> ,N=i,l5> 

?ITe!:(6,Fl0)  ( 7(4, M)  ,N=1,15> 

RITE(6,5U)  (7(6, N)  ,N=1,15) 

5TTC  1 1.  5 191  I 7 I Ml  . M = i . i R\ 


WRITE(6,*13)  (7(7, N)  ,M=1,  15) 

CONTINUr 

Ie(IlN  .E".  1)  RETURN 


icoun=o 


(3- 


51  icoun=ico'in*i 

IF(MM-TCHS(  TG)  ) ^501,3500,3601 

3600  IP=I»+1 
TCOMN=TPRM(IP) 

3601  CONTINUE 
RETURN 
END 

C ///////////////////////✓////// ///////////////////////////////////////// 
C 

SUBROUTINE  REYSTR(<ON,T3, X,TREF,XNUE,X3E,S,ITCNT1, 

1RCON, ALPHA.ITRAN) 

c 

c ///////////////////////////////////////// // ( n f r f t r r f rr  u t r / r r r 1 1 1 f / rr / 

COMMON  3,  PR,  RFY,  XMINr,  OMEGA,  00,  TW,  P10,  T10,  RIO,  VIS10,  TE, 

1 PEfREjJEjVIPINF^U^oSt^OYWjSItERRORiTCfTA^rOGEjIENDl.INTACT, 

2 PRT,XX<, 1TRX,XLAM,VA  RPRT  , XI  NT  ER  ,SE°0 , 1 CHS  ( 3)  , I PEN  (3), *0(20  0)  , 

3 EN(2  0 0)  ,EP(200)  ,ETO(200>  ,ETN(  2 00>  , ET 5 ( 20 0)  , FO  ( ? 00 ) , FN ( 20  0)  , J2D A, 

4 rP(  20  0)  , T\|  ( 20  0 ) , TO ( 2 0 0 ) , XNN  (20  0 ) , VM  ( 20  0)  , VO  (70  0 ) , VP  (20  0 ) ,TP  ( 20  0)  , 
501(?0  0),'>?(200),03(200),T?JP(200>,T7JO(200),T2JN(200> 

TTR=  ( T A+ 1 12.  )/(TA*TRFP4-il2.) 

CO=T  0 ( 1) 

00=FPM)  =XNN(1)  =TPI  = 6LT=n. 

C SHEAR  STRESS  AT  THE  WALL  AS  THE  SCALING  FUNCTION 

Yll=  ( (2.  4-XXO  * (1,  ♦XXK*XXk'**2)+i.+XX<>  /(  (l.+XXO  Ml,  »-XXK  + XXK**2)  ) 

Y 12= ( 1 . ♦XX<*XX<**2)/X<<**2 
Y13=(l,fXXKfXXK**2>/ (XX<**3* (l.+XXK) ) 

Y 14=1,  /(XXK*V3*  (l.+XXK+XXK**2)  > 

FETW=(-Yll*-D( 1) +Y12*rP(2) -Y13*FP(3) M 14*Ft>(4) ) 70YW 
FFTW=ARC(cErW) 

XLMlW=((l,+rR)*SORT(TP(l) )/(TP(l)fTR) ) 


»I2=XL  *1  M*'rrTW 

SHEAR  STPESS  AT  THE  WALL  AS  THE  SC  ALIMS  FMNCTIOM 
0 0 1 N=  ? , < O M 
OY=OYW*VXK** (M-2) 

XLM1=  < <1.  +TR)  *SQRT(TP(N)  ) /<TP(N)  +TR)  ) 

C=TD(M) 


TDI=TPI+.  9*0Y*  (CO+C) 

cn=n 

XNN(N)  =roST3l*?ORT(2.  *X)  /(RE*U£) 

I p ( J’O  A.  ME.  0 ) XMN  (M)  =XM'I  (N)/RCON 

TF(TT9AM,  ME.  0)  XMN(M)  = (:>rnM/COS(ALPHA)>*(-l.+S9PT(l,+ 

12.*S0pt{?. *X) *E°S*COS  < ALPHA) *TPI/(RE‘JE*RCON**2) ) > 

TP(nLT.r’T.  0.  ) GO  TO  2 

Ic(pp(M> . GE. 0. 9Q5)  3LT  = XMN (N)-  <FP(M)-.995)*<XMM<N)-XNN(N-1)) 

1 /(FP(H)-po(N-I  ) ) 

IPdTRAM.  "0.  i)  GO  TO  1 2 

00=00+  ((1.  -PD  (M)  )*TP(V)  (l.-FP  (N-l)  ) *TP  (N-l)  ) *0Y/2. 

12  CONTIMIP 

I F ( IT p AM, ME.  0)  00=00+ ( (l.-FP(M) >*TP(M> /SORT ( T2 JP(M)  ) +(l.-FP(s- 

11))*TP(M-1)/S0RT(T?JP(N-1 ) ) ) *0Y/2. 
oil=SOPT(->,  *X«RPY/(TREF»*i,5*TTR))  * T°T*  * 2/ (XNUE  *TP  ( M ) **  3 ) 

I c ( J20 A. ME.  0 ) PI1=PI1YROON 

Ic( TTRAM. ME, 0)  PI1= ( (RE*UE) **2> * (XMM(M) **?) »PS0N*?9RT<T2 JP(M) ) 

1*  <REY**1.  G)  / ( (TP(N)  **3) “SORT  (2  . * X)  * ( TRE  F**  2.  25 ) * ( TT  R 
2**1.  9)  *XM'IE) 

OY=OYW*XV<** (M-l) 

0 YM1=0 Y/XXK 

Y8=2./  (fOYMl/OY)*  (1.  +OYM1/OY)  ) 

Y9=2.  / (1.  +OY/OYH1) 

Y10=l. -0Y/0YM1 

PT2=XLM1*EP(M)*  AOS(Y9»FP(N  + l)/2.-Yn*-P(N)-Y8*rP(M-l)/2.  ) /DY 

CEBIOE-SMTTH-MONSIMKIS  "00Y  VISCOSITY  MOOEL 
Y°LUS=Sr'DT(DIl*PI2)/  (?S.  *XLH1) 

IPCYPHI9.GT.5n.)  YPLUS  = 'P. 

EP(N)=«1G*PI1* ( 1. -EXP C-v PLUS) ) **2*A0S(Y9* 

1 FP(N  + l)/2.-Y10*FP(M)-YA*FP(N-l)/2.)  /(0Y*XLH1) 
CE0ICE-9HTTH-H0NSIMKIS  rO0Y  VISCOSITY  MODEL 

TRUNCA  Tr  THE  IMMER  RESIOM  CALCULATIDM 
IFCEO(N)  ,LE.  E°(N-1) ) EP  ( M)  =E  P ( N-l) 

TRUNCATE  THE  IMMER  RESIOM  CALCULATION 

CONTINUE 
00  3 N=1  ,KOM 

XLM1=  ( (1  . +T  R ) * S ORT  ( T°(  M)  ) / ( T P ( N)  +T  R)  ) 

001=. 0199*9 ORT ( 2.*X*REY' (TREF**1.5*TTR) ) *DDJ (XNUE* XlMI*T PCN) **2) 

I e ( J20 A,  ME,  0 ) OnirOOlYROOM 

IFCDOl.LE.EP(M) ) EP ( M ) =001 

ARGE=-.M?*<  (S-9TPX)  / XLAH)  **2 

IFCARGE.LT, -57E. 83)  GO  TO  9 

E°(N) =FD (M) * ( 1 . -EX® ( A ?Gr ) ) 

9 CONTINUr 

IFCXINTPR.EO.n.)  EP(M)=1.+EP(N) 


4 


fl 

t: 

S 

I 

I 


1 


IC(XINTF=>.F0.  1.  ) FP(M)=1  . ♦ (1. 75/(1. ♦5.5*(XMN('D/9LD**6)*1.)* 

1 FP(N)  /?.75 

3 ETP(N> =1. +RR* (EP(N)-l. > /PRT 
C 

RETURN 

END 

c////// ////////////////////////////////////// /////////////////////////// 
c 

SUPPOUTT^  MATE "NT  (Xl,Y?,X3,Yl,Y2,Y3,411,At?,A13,A21,A22,A23, 

S A31,A3*, A33,LC,LN,L0> 

c 

c ///////////////////////////////////////////////////////////////// ////// 

c 

C THIS  StnPRjrrNE  SOLVES  SHE  THREE  SIMULTANEOUS  3 AMR  MATRIX 

C EOUATTCMS 

c 

C Atl*Xl  «■  A1?*Y2  + A 1 3 * X 7 = Yt 

C A 71*X1  f A2?*X?  ♦ A23*X7  r Y2 

C A31*Xi  f A32*X2  + A33 *X3  = Y3 

C 

C POR  XI,  X2,  AMO  X 3 

C 

C A TJ  ARE  9 SAND  MATRICES  OF  LENGTH  LO,  WORKING  LENGTH  LN, 

C AND  WI 0TH  LC 

C (THESE  MATRICES  ARE  ASSUMED  TO  9E  CORNER  ADJUSTED,  I . E.  THE 

C CORNER  ELEMENTS  ARE  STORED  IN  (1,1)  AND  (LN,LG),  ETC.) 

C 

C XI  AND  YT  ARE  VECTORS  OP  LENGTH  LQ  AND  WORKING  LENGTH  LN 

C 

C 

DIMENSION 

* Xl(LD) ,X?(LD)  , X3(LD) , Yi ( LO) , Y2 (LQ) , Y3 ( LQ) , 

* A11(LO,LO  , A12(LQ,LC)  ,A13(LQ,LC), 

* A21 (L 0, L 0)  ,A?2(L0,LC)  ,A23(L0,LC), 

* 431 (LO,LC) , A 72 (LO,LC) , A33(LQ,LC) 

C 

C INITIALISATION 

C 

c 

L °=LN+ 1 
L = (LC-  1)  /? 

LM=LN-L-1 


IF(LC.GE.LN)  L=LN 
00  3 1=1  , LN 
X1(I) =Y1 (I) 


X2(I)=Y?  (I) 

X3(I)=Y3(T) 

3 CONTINUE 
C 

C DOWNWARD  GAUSSIAN  ELIMINATION  WITH  =>IVOTING 

C 

c 

00  401  K=1,LN 
IF(L.EO.LM)  L = LN 


IF(l.LT.LM)  t=L  ♦ 1 


U=AnS(All(K,l>> 

I ~K 
H=1 

no  113  J=<»  L 
IPU.'-O.K)  00  TO  111 
v=Ans(Aii(j,n ) 
IF(V.L'.H)  GO  TO  111 
U=V 
M = 1 
T=J 

111  V=A9S(A?1  <j,m 
TP(V.LP.U)  30  TO  112 
U=V 

M=2 

I=J 

112  V=A9S(An(J,l>> 

I F ( V • L “•  'J  ) GO  TO  113 

(J=V 

*1=3 

T=J 

113  OOMTINU* 

IP(I.Sn,<)  30  TO  115 
IP(M.MP.t)  30  TO  116 
00  114  J=  t»  L 3 
U=A11(K, J) 

A11(K,  J)=A11(I,  J) 

All  (I*  =') 

U=A1?(K, J) 

A12(K,  J)  = U2(I,  J) 

A 12(1 1 J)  ='J 
U=A13(K, J) 

A13(K,  J)=A13(I,  J) 
At3(I,  J)='l 

114  CONTINU- 
U=X1(K) 

X 1 (K) =X1 (I) 

X1(I)=M 
GO  TO  120 

115  IF(M.PO.l)  GO  TO  120 

116  IF(M.N5.?>  30  TO  113 
00  117  J*1,L3 
U=Alt(K,  J) 

A 11 (Kt  J) = A71 ( T»  J) 

A 21  (I,  J)  = 'J 
U=A12(K,  J) 

A12(K,  J>  = A2?(I,  J) 
A22(I,J)=U 
U=A13(K,  J) 

A 13( K*  J) = A23  (I  t J) 
A23(I,  J)=H 

117  CONTINU” 

U=X1(K) 


X1(K)=X?(I) 

x?m=’j 

GO  TO  l’O 

11*  00  110  J=1,LC 

U=A11  C K.  J> 

A11(K,  J>  = 131  (I,  J) 

ATI  (It  J)  ='l 
U = A12(K,  J) 

A 12( K,  J)  = AT?  ( I , J ) 

A 32(1,  J)  ='l 
U = A13 ( K, J) 

A13(K,  J)  = 133(I,  J) 

A33(I,J)='I 
110  CONTI  Nil; 

!J=X1(K> 

XI  (K>  = X3  ( T) 

X 3 ( I ) -U 
120  CONTINUE 
C 

00  128  1=  <»  L 
TF(I,FO.K>  30  TO  123 
U=A11  (I*  1)/A11  ( K,  1) 

00  122  J = it L 0 

IFtJ.NM)  A13(I,  J-1)=M1  (I*  J)  -All  (K,  J)  *U 
A11(I,  J>  = 112(1,  J)-A1?(K,J)*U 
A12(I,J>  = A13  (I,  J)-A13(K,  J>*U 

122  CONTINUE 

A 13(  I » L'*)  -=0# 

xim=yim-xi  m*u 

123  CONTINO" 

U=A21(T,1)/A11(K,1) 

00  120  J=t,L3 

IF(J.Nc.  1)  A23(T, J-l>  = A',1(I,J)-A11(K»J)*U 
A 21(1, J)  =422(1, J)-A12(K, J)*U 
A22(I,J)  = A?3(T»J)-A13(K,J)*U 
125  CONTINUr 

A 73(1 , Lr)  =0. 

X2(T>  = X?(T)-X1(K)*U 
U=A31(I,1)/A11(K,1) 

00  127  J= 1, lC 

I F ( J,  Nr«  1 ) A33(T,  J-l)  =A’l  (I,  Jl-All  (<,  J)*U 
A 31 (I,  J)  = 132(1, J)-A12(K, J)*U 
A32(I,  J)  = 133(1,  J)-A13(K,J)*U 

127  CONTI MM* 

A 33 (I, LO)  =0. 

X3(I)=X*(I)  -XI  (K)*U 

128  CONTINtr 

c 

U=AOS ( A21 (<, 1) ) 

1 =K 
M=2 

00  213  J=K,l 
IF(J,FO,K)  GO  TO  212 
V=AOG( All (J, 1) ) 


' 


IP(V.LF.'J)  GO  TO  211 

U=V 

M=1 

T=J 

211  V=AOS(A2i (J,l> ) 
IP(V.LE.'J>  GO  TO  21? 

u=v 

H=2 
I = J 

212  V=AO<; (ft71  ( J* 1 ) ) 
IP<V.Lr.IJ)  GO  TO  213 

u=v 

M=3 

I=J 

213  CONTINfr 
TF(T.FO,<>  30  TO  215 
I F CM.  M "#  2)  GO  TO  216 

no  ?i4  j=i,lo 

U=A2t<«,  J) 

A21(K, J>=121 (I, J) 
A2KI,  J)='J 
U=A2?(  K,  J> 

A?2(K,  J)  = 42?  ( I » J) 

A 22(1, J)  = U 
U=A23(K, J) 

A23CK,  J)=A?3(T,  J) 

A 23(1,  J>='J 

214  rONTTMir 
U=X?C<) 

X2(K)=X2(T> 

X 2 ( T ) =0 

GO  TO  2?n 

215  IF(M,N“,3)  30  TO  220 

216  IF(M.N'.l)  GO  TO  218 
00  217  J = l,  LG 
ll=A21  ( K,  J) 

A 21  (<,  J)  = Ml  (I,  J) 
A1KT,  J)  =>» 
tf=A2?(K,  J> 

A22(K,  J>=A12(I,  J> 

A 12(1 , J)  =0 
U=A?3(K,  J> 
A?3(K,J>=M3(T,J) 

A 13(T, J>  *M 

217  CONTINUr 
U=X2 (K) 

X?(K)  = X1(T> 

xim=u 

GO  to  22  0 

218  00  219  J=1,LC 
tl=A21 (K, J) 

A21(K,  J)*A31  (I,  J) 

A 31(1,  J)  ='l 
U=A22 (<, J) 
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■ 


A??(K,  J>  = A3?  <T,  J) 

A3?(T,  J)  ='» 

U=A?3(K,  J) 

A’3<«,  J)  =AT3(I,  J) 

A 33(1,  J)  = 0 
?m  C0NTiMHr 
U=X2(K> 

X?(K)=X3  (T) 

X3(I)  = U 
220  CONTINUE 

00  225  I =X,  L 
Te(T.EO,<)  50  TO  223 
U=All(r,t)/A2t(K,l) 

00  22?  J=1,.0 

I c ( J»  NET.  1 ) A13CI,  J-l)=Ail  (I,  J)  -A21(K,  J)  *U 
A1KT,  J)  = A12(I,  J)-A22(K,  J)*U 
A12(I,J)  = H3(I,J)-A23(K,J)*U 
22?  CONTI  HO: 

A13(T,L2>=0, 

X1CI)=X1 (I)-X?(K)*U 
ll=A21<I,l)/A?l(K,l> 

00  2?5  J=l,LC 

IP(J.NC.  A?3(T,  J-l)  =AM  (I,  J)  -A21  (<,  J)  *U 
A?1(T,J)  = 1?2(I»J)-A?2(K,J)*U 
A‘*2{I»J)=A2  3(I»J)-A?3(K,J)  *U 
225  CONTI  NIJP 

A'»3(I,LO)=0, 

X?(I)=X?(T>-X2(K>*U 
223  CONTINUE 

U-=A31  <I,1>/A21  (K,l) 

no  227  J=l,l5 

IP(J.Nc.  1)  A33(I,  J-1)=A71(I,  J)-A21(K,  J)  *U 
A3KI, J)=A3?(I, J)-A?2(K, J)*U 
A 32(1, J)  = A33(I, J)-A23(K, J)*U 
227  CONTI NOr 

A 33 ( I , LO) =0, 

X3(I)=X-*(I)-X2(K)*U 
225  CONTINUE 

IF(K.PO.LN)  50  TO  401 
U=ARS (A7t (K, 1) ) 

I=K 
M=3 

JL=K*1 

00  313  J=JL,L 
V=AR5( All  (J, 1) ) 

IP(V.LC.U)  GO  TO  311 

u=v 

H=1 

r=j 

311  V*A3S(A?1<J,1>> 

IF(V.L5.0)  50  TO  312 
U«¥ 


.^1 
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I=J 

11?  V = A«S(  AM  (J,  1>  > 

IF(V.LE.U)  30  TO  313 
ll=V 
M = 3 
I=J 

313  CONTINUE 
TF(I,rO,<)  GO  TO  320 
IF(M.M".3)  30  TO  316 

no  314  j=i,  - : 

ll=A31  (K,  J> 

A31K,  J)=A31(T,  J) 

A 31(1,  J)='» 

U=A?2(K, J) 

A3?(K, J)  = A32(I, J) 
A32(T,  J>='J 
U=A33(X,  J) 

A33(K,J'  = A33(I,J) 

A 33(1,  J)  ='» 

314  CONTINUE 
IJ=X3(K> 

X3(K)-X3  ( T) 

X3(I)='J 

GO  TO  3?1 

316  IFtH.MF.il  30  TO  313 
00  31’  J= 1, L 3 
U=A31 ( K,  J) 

A31(K,  JI*U1(I«  J> 

All  (I,  J)  = ') 

IJ=A32(K,  J> 

A32(K,  J)  = A1?(I»  J> 

A 12(1,  J)='J 
U=A33(K, J) 

A33(K,  J)=U3(I,  J) 
A13CI,  J)  = U 

31’  CONTINUE 
U=X3(K> 

X3(K)=X1(T) 

Xi(I)=U 
GO  TO  3?0 

31ft  00  319  J=1,L3 
U=A31  ( K,  J) 

A31(K,  J)  = 421  (I  , J) 
A?1(T,  J)  =U 
tl=A32(K,  J> 

A32CK,  J)  = 422  (I,  J) 
A?2(I,  J>=U 
U=A33(K, J) 

A33(K,J)  = 423(I,J) 

A 23(1,  J»  =U 

319  CONTINUE 
U=X3(K) 

X3(K)  =X?(I> 


o o o <r> 


x?m=u 

320  rONTINU- 

c 

Tl=««-1 

no  3?a  t=il,l 
•J=All(T,l)/A3i(K,l) 
no  32?  J=l,L3 

IPCJ.Nr.  1)  a 1 3(  I,  J-l)  = A11  (I,  J)-A31f  K,  J»  *u 
A 11  Cl,  J>  = At?(T,  J)-A3?(K,  J)*'J 
A 12 (I, J)=A13CI, J>-A33(<, J>*U 
32?  CONTINUE 

A13CI ,L?) =0. 

Xl(I>  = X1  (T)-X3(K)*U 
U=A21(T,1)/A31(K,1) 

DO  3?*!  J=l,L0 

IF(J.NE.l)  A 23(1,  J-l)  =A’l  (I,  J>-A31<K,  J)  *U 
A 21 (I, J>=A22 (T, J)-A32(K, J)*U 
A22CI, J>=A23(I, JJ-A33CK, J)*U 
325  CONTI  MiF 

A 23  ( I , L"*)  =0. 

X?(I)=X?(T>-X3<K>*U 
U-A31 (I,1)/A?1(K,1) 

00  327  J=1,LC 

TPCJ.NE,  t>  A 33  (T,  J-l)  =A*1  (I,  J»  -A  31  <K,  J)  *»i 
A 31 (I, J)=A32(I, JJ-A32CK, J)*U 
A3?(T, J)=A33(I,  J)-A33CK, J)*U 

327  CONTINUE 
AISCTtL'M^O. 

X3(I)=X-»(I»-X3(K)»U 

328  CONTINU- 

c 

401  OONTTNU- 

UPWAPD  0AMS5IAN  ELIMINATION 


L =1 

00  507  X=1,LN 
I=LP-K 
C 

U=X3 ( I > 

IPCT.PO.LN)  30  TO  502 

00  501  J=?,L 

IJ=H-J 

501  (J=U-A3?(T  , J-i) *X1  (IJ-1)-  A 33(1,  J-l)  *<?(  I J-l>  -A3t  ( I,  J)  *X3(  IJ-1> 

IPCL.O-.LO)  U=U- A 32  ( I , LC)  *X1  ( I +LC)  - A 33  ( I , LSI  *X?  <I*L  Z ) 

50?  X7(I>=U/A31(I,1) 

C 

U=X2(I)-A?2(I,1)*X3(I) 

IPCI.EO.LN)  GO  TO  504 
00  503  J=?,L 
IJ*I+J 

503  U=U-A2  3(T,J-1)*X1(IJ-1)-A21(I,J)*X?(IJ-1)-A22(I, J>  * <3C I J-l> 
IFCL.GE.LO)  UslI-A23(I,L0)*Xl(I+LC) 
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ooor»oo  o r> 


504  X’(I)=M/A21 (1,1 ) 

C 

u=xici)-ai2(I,i>*x2(I)-ai3(i,i)*x3(I> 

IPd.^n.L'O  GO  TO  505 
no  505  J = ?,L 
I J=I *J 

505  U=»J-A11(T,J)  *X1  (IJ-1) -A12 (I, J)»X2(IJ-1) -A13(I,  J)*X3(IJ-1) 

505  X1(I)=M/A11(T,1> 

TP(L.LT.LC)  L =L  *1 
C 

507  CONTINUE 

c 

RETURN 

end 

//////////////////////////////////////////////////////////////// /f///// 
SUBROUTINE  CONST  ( ICON,  I5EF) 

/////////////////////////////////////////////////////  ru  nr  ft  / rf  n ///// 
ICON=0  ( A XI SYMMETRIC  30TY  RAOIUS  IS  COMPUTEO  FROM  SJRFACE 
LOCATIONS) 

1  (RADIUS  POLYNOMIALS  ARE  IN°UT  TO  THE  PRDSRAM  BY  REAO 
STATEMENTS  LATE5  TN  THIS  SUBROUTINE) 

COMMON/'ONST/  R ( 30)  , 7 ( T 0)  , XP  ( 3 0)  ,A  ( 3 0)  , 9 ( 3 0)  , C ( 70)  , 0 (30 ) , EC  30  > , 

IF  (30)  ,C=>  (30)  ,OP  (SO)  ,INUM 

COMMON  5,  or,  °EY , XMTNr,  OMEGA,  90,  TH,  310,  T10,  R10,  VIS13,  TE, 

1 PE,RE,U- ,VISINF,SU,"oS, DS,OYW, SI, ERROR, TC,T A, I EDGE, IEND1, INTACT, 

2 PRT,XXX,0TRX,XLAM,\/ARP5T,XTNTER,SEP0,ICHS(8),rPRNn>,E0  (20  0)  , 

3 EN(20  0),EO(2n0),ET0(?0»,),ETN(?0  0),ETO(20  0),F0(?0T),FN(?n0),J20A, 

4 CP(200) ,TN(200) ,TO(200) , XNN  ( 2 0 0 ) , VN ( 20 0) , VO ( 20 0 ) ,V*(230),T®(200), 
501(200), 02(200)  ,03(200)  ,T?JP(200),T2JO(200)  ,T2JN(200) 

I F (ICON#  rD,  1 ) GO  TO  25 
WRITE  (5,  2 41) 

241  FORMAT (5X,“R/L  7/L") 

OO  32  IT  si, INUM 

WRITE (5, 9 01 ) »(IT)tT(IT) 

801  FORMAT (?ri5. 9) 

32  CONTINUE 
<END=I NJM-2 
N=1 

00  19  K=l,KENO 
NEN0=N*2 
00  12  M=N,NEND 
0 (M) =1. 

E (M) =2  (M) 

12  F(M)=7(M)**2 

OET  =0  (N)  *E(N«-1)  *F(N»?)  *E(N)  *F  (N+l)  *D(N+2>  fO(NH)  *£(N*2) 
1*F<n)-F(N)*E(NM)*D(N«-2)-E(N>*0(N*1>  *F( N+2) -D( N)  * 

2F (N+l) *E (N+2) 

00  14  M=N,NENO 
lh  P(M)=R(M) 

0ETA=0(N)  *E  ( N+l  ) * F(  N*2)  *F(N)  *F  ( N*l»  • 0 ( M «-2)  *0  (N*  1>  *E  ( N*2 ) 
1*F(N)-F(N)*E (N*1)*0(N*2) -E(N)*n(N«-l>  *r(N*2)-D(N)  * 

2F  (N+l ) *ECH2) 
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DO  1*  M=N»N"MD 
D (H)  =1,0 
15  F(M)=o<M) 

0FTP=0  ('I  > *E(N4-1)*F(N4-?)  «-F  (N)*F  (N  + l>  *0(  N«-2)  4-0  < 1 ) *"(N«-2> 

1*F(N)-F(M>**  (N+t)*D(N*?>  - "(N)*0(Nf  1)  *F(n+2I-0(M)  * 

2F  (N+l)  *E  (N*2) 

00  19  «=N,NENO 
E CM)  =7  (M) 

19  F(M)=P.  (M) 

OETC=D(N)*F  <N+l)*FC'l  + ?)+rm*r  (N+i ) * 0 CM  f?)  + 9 ( M *1 ) *E ( N+2) 
i*F(N)-F(N>*z  (N  + l)*9(Nf2)-E(N>*0(N4-l>  *-  { N + 2) -0  ( M ) * 

?F (N+l) *E(N*2> 

A (t<)=OrTa  'OET 
0 (K)=OrTq/orr 
c (K)  =rtCT':  /OET 
N=N*1 

10  CONTTMDr 
7CON=0.0 
RCON=0.P 
0=0.  0 
N=1 
KN=1 
07=00 

ZZP.FFsZCTNU'C)  +0Z 
22  CONTINUE 

IF(A0S(7CON-7(KN> ).LE.07/2.)  XP(KN)=S 
IFCKN.EO.TNUM)  GO  TO  25 
TP(A00(7OON-?(KN) I. IE. 07/2, ) KN=KN+1 

25  CONTIMUr 

IFCN.EO.  <-MD)  GO  TO  24 
ZPFF=Z  (N*  1) 

ifczcon.gz,zr:f>  n=n*i 

24  0P07=O(M) *?. 0*C(N)*7CON 
ALPHA=ATAN(OROZ) 
nP=nS*0TM(ALPMH) 

OZ=DS*  000 (ALPHA) 

RCON=P.COM*OR 

S=S*O0 

ZCON=ZrOM^OZ 

IFCZCON.LT. Z2PEF)  GO  TO  22 
no  57  KL=1,IMUM 
0 (KL)  = 0. 

F CKL)  = 0. 

83  CONTINUE 

26  CONTINUE 

I F (ICON.  FO.  0 ) GO  TO  41 

J=1 

K=i 

RCOM=0. 

ZCON=0. 

S=0. 

K N=1 
D7=0S 

RFA0C5  f 1 000)  IREF 
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liTon  format (IT?) 

1001  FORMA”-  ('fi5, 9) 

100?  FOPMftT  (1^15. 9) 

00  ??  1=1  * IRE F 

RFAn(5,10  01)  A(T),n(I)t'*<I),D(I),ECI> 

RFA0(5»inn?>  F(i) 

27  CONTimi* 

fCIREf>=f(IR-p>  +02 

36  IPtAOSC^OON-ZKNn.LF.OT/?.)  XP(KN>=S 
IF(KN.  "•’.TMJH)  T,0  TO  ?9 
IF(A90(700N-7(KM)  l.LE.OT/?.)  KN=KN*1 
29  nPnz=9(J)  f?.*C<  J)  *700N4-3.*0(  J)  *ZC0N**?+4,*E< J> *700N**3 
ALPHA= Ar AMtOROZ  ) 

0‘,=0S*COP  ( AL  pHA ) 

s=s*os 

7CON=700N*-r)7 

pCON=A (J)  +3(J) +7C0N  + C (J) *7C0N*»2*D< J > * ? CON** 3f E ( J) *ZC0M**4 
IF(700N.OE*- ( J)  ) J=J*1 
IF(7C0M.LT,r(i^CF))  GO  ?0  36 
WPIT E(  6,  ? '*?) 

242  FORMAT  (inx, "C®  7/1  **) 

00  3 < 9=1  ,INU*-« 

H?ITF(6,’U?)  0 ° ( K9)  » Z ( <9 ) 

7112  FORMAT f’F15. 9) 

3 CONTINUE 
41  CONTTNU' 

RETURN 

fno 


105 


Vita 


Charles  R.  Blake  was  born  October  22,  1946  in  Meridian 
Mississippi.  He  received  his  primary  education  in  Hickory, 
Mississippi  and  graduated  as  valedictorian  from  Hickory 
High  School  in  1964.  He  received  his  undergraduate  education 
at  East  Central  Junior  College  and  Mississippi  State  Univer- 
sity. He  was  awarded  a Bachelor  of  Science  Degree  in  Aero- 
space Engineering  from  Mississippi  State  University  in 
May  1969. 

Upon  graduation  from  college,  Capt  Blake  was  employed 
as  an  aerospace  engineer  for  Brown  Engineering  Company, 
Huntsville,  Alabama  until  entering  the  USAF  in  August  1969. 
He  received  a commission  in  November  1 969  from  the  USAF 
Officer  Training  School  and  was  then  assigned  to  Under- 
graduate Pilot  Training  at  Moody  AFB,  Georgia.  He  received 
his  wings  in  December  1970. 

Capt  Blake  then  served  as  a C-130  pilot  and  flight 
instructor  at  Little  Rock  AFB,  Arkansas,  CCKAB,  Republic 
of  China,  and  Clark  AB,  Phillipines.  He  was  then  serving 
a second  tour  at  Little  Rock  AFB  immediately  prior  to 
entering  the  School  of  Engineering,  Air  Force  Institute  of 
Technology  in  May  1975. 

Permanent  address:  Route  1, 


Chunky,  Miss.,  39323 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  fWhm  Dele  Entered) 

REPORT  DOCUMENTATION  PAGE 

t.  REPORT  NUMBER  / |2.  GOVT  ACCESSION  NO. 

GA/MC/76D-U  / 

«.  TITLE  (end  Subtitle) 

NUMERICAL  SOLUTION  OF  THE  COMPRESSIBLE 
BOUNDARY  LAYER  EQUATIONS  OVER^^^ 
AXISYMM2TRIC  SURFACES 


i press: 


17.  AUTHORS 


Charles  R.  Blake 
Capt  USAF 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Air  Force  Institute  of  Technology  (AFIT-EN 
Wright  Patterson  AFB,  Ohio  ^5^33  s 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 


4.  MONITORING  AGENCY  NAME  & ADDRESSf//  dltlerent  trom  Controlling  Ol/ice) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 
3.  RECIPIENT'S  CATALOG  NUMBER 


3.  TYPE  OF  REPORT  6 P*R*B©-CO  VERED 

MS  Thesis 

6.  PERFORMING  ORG.  REPORT  NUMBER 
8.  CONTRACT  OR  GRANT  NUMBERf*.) 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  8 WORK  UNIT  NUMBERS 


12.  REPORT  DATE  * 

Dec,  1976 

13.  NUMBER  OF  PAGES 
106 

IS.  SECURITY  CLASS,  (ol  thle  report) 


Unclassified 


ISa.  OECLASSI  FICATION/ DOWNGRADING 
SCHEDULE 


118.  DISTRIBUTION  STATEMENT  (ol  thle  Report) 


Approved  for  public  release;  distribution  unlimited 


J *7.  DISTRIBUTION  STATEMENT  (of  the  ebetrect  entered  In  Block  20,  i(  different  from  Report) 


I U.  SUPPLEMENTARY  NOTES 


IS.  KEY  WORDS  ( Continue  on 

Boundary  Layer 
Compressible 
\ Axi symmetric 


^pp^vedAf^^uMic  release;  IAW  AFR  190-17 

Jei(r£l  F.  GuessT^ Captain,  USAF 
Director  of  Information 


eery  end  identify  by  block  number ) 


20.  ABSTRACT  ( Continue  on  reveree  elde  If  neceeeery  end  Identify  by  block  number) 

A numerical  solution  of  the  compressible  boundary  layer  equations 
was  developed  for  flows  over  either  two-dimensional  or  axisymmet- 
ric  surfaces.  The  solution  method  is  capable  of  solving  for  bound- 
ary layer  parameters  in  either  laminar  or  turbulent  flows.  In  the 
case  of  turbulent  flow,  closure  is  achieved  by  use  of  a two-layered 
eddy  viscosity  model.  The  boundary  layer  equations  are  solved  by  a 
numerical  marching  procedure.  A Mangier- Levy-Lees  transformation  of 
independent  variables  is  used  to  improve  the  efficiency  of  the\ 


DO  1 JAN*7J  1473  EDITION  OF  1 NOV  8S  IS  OBSOLETE 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dm r.  » ■ td) 


